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Abstract 


Stagnation  point  flow  and  heat  transfer  in  the  presence  of  free-stream  tur¬ 
bulence  is  investigated  through  both  numerical  simulation  and  theoretical 
analysis.  Large  eddy  simulations  (LES),  using  fourth  order  finite  difference 
in  curvilinear  coordinates  and  an  efficient  dual  time  linearized  sub-iteration 
scheme,  are  performed  to  investigate  free-stream  turbulence  impingement 
upon  an  elliptical  leading  edge  and  the  resulting  heat  transfer  enhancement. 
A  new  blending  procedure  is  developed  through  which  independent,  statisti¬ 
cally  identical  realizations  of  homogeneous  isotropic  turbulence  are  combined 
to  provide  realistic  representations  of  free-stream  turbulence. 

Results  for  different  free-stream  turbulent  intensity,  length  scale,  and 
Mach  number  are  reported.  Strong  anisotropy  of  the  turbulence  due  to  the 
mean  flow  strain  is  observed  as  the  stagnation  point  is  approached.  The 
Reynolds  stress  statistics  and  budgets  are  obtained  and  presented.  These  re¬ 
sults  are  expected  to  provide  unique  data  for  turbulence  modelling  of  strain 
dominated  flows.  The  numerical  results  show  good  agreement  with  the  ex¬ 
perimental  measurements  on  elevated  heat  transfer  coefficient.  It  also  reveals 
that  small  scale,  intense  vortical  flow  structures  generated  at  the  leading  edge 
by  vortex  stretching  induces  significant  changes  in  local  thermal  boundary 
layer,  causing  the  observed  heat  transfer  enhancement. 

In  the  theoretical  study,  the  distortion  of  three  dimensional  unsteady  dis¬ 
turbance  in  an  incompressible  Hiemenz  boundary  layer  and  its  effect  on  the 
wall  heat  transfer  is  analyzed  based  on  linear  vortex  dynamics.  An  asymp¬ 
totic  solution  for  disturbance  vorticity  is  obtained  with  explicit  dependence 
on  the  disturbance  length  scale  and  frequency.  It  is  shown  that  the  vortic¬ 
ity  amplification,  and  hence  the  heat  transfer  enhancement,  increases  with 
decreasing  length  scale  and  the  maximum  value  is  found  around  five  times 
the  boundary  layer  thickness.  The  unsteadiness  of  the  disturbance  reduces 
the  disturbance  amplification  but  is  of  second  order  at  low  frequency.  By  ex¬ 
tending  the  analysis  to  free-stream  turbulence,  a  new  scaling  correlation  for 
the  relative  heat  transfer  enhancement  is  derived  which  incorporates  turbu¬ 
lence  intensity,  integral  length  scale  and  mean  flow  Reynolds  number.  This 
correlation  is  shown  to  reasonably  collapse  the  recent  experimental  data. 


Keywords;  Stagnation  point  flow,  free-stream  turbulence,  turbine  blade 
heat  transfer,  large  eddy  simulation,  vorticity  dynamics,  turbulence  mod¬ 
elling. 
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Chapter  1 


Introduction 


1.1  Motivation  and  Background 

Stagnation  point  flow,  where  fluid  approaching  a  solid  surface  divides  into  di¬ 
verging  streams,  is  a  ubiquitous  and  important  flow  phenomenon.  Hiemenz’s 
pioneering  work  at  the  beginning  of  last  century  showed  that  the  incom¬ 
pressible  stagnation  point  flow  is  one  of  those  very  few  fundamental  types 
that  admit  exact  solutions  to  the  Navier-Stokes  equations.  However,  this 
seemingly  simple  flow,  when  disturbed,  exhibits  surprisingly  rich  phenomena 
and  intriguing  dynamics.  This  can  be,  at  least  partially,  attributed  to  the 
fact  that  in  the  stagnation  region  the  upstream  disturbance  typically  reaches 
its  maximum  by  the  effect  of  mean  flow  strain  while  the  mean  flow  velocity 
itself,  on  the  other  hand,  approaches  zero  —  perturbed  stagnation  point  flow 
is  inherently  nonlinear.  This  peculiar  feature  contributes  to  many  difficulties 
in  our  understanding  of  the  nature  of  the  perturbed  stagnation  point  flow. 
The  correct  account  for  the  linear  instability,  for  instance,  did  not  appear 
until  late  1970’s. 

In  addition  to  the  substantial  theoretical  interests,  the  practical  reason 
why  stagnation  point  flow  continue  to  remain  as  an  active  research  subject 
after  more  than  a  century’s  study  lies  in  its  ubiquitousness:  almost  all  the 
interactions  between  fluid  flow  and  solid  structures  involve  some  kinds  of 
stagnation  points  or  lines.  Many  flow  and  heat  transfer  problems  encoun¬ 
tered  in  various  important  engineering  applications  are  of  stagnation  point 
nature.  Gas  turbine  cooling  and  micro  electronics  design  are  among  the  most 
noticeable  examples. 
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As  one  of  the  major  energy  conversion  devices  of  our  times,  gas  turbine 
is  used  both  for  aircraft  propulsion  and  in  land-based  power  generation.  Up 
to  now,  it  is  widely  recognized  that  the  key  technology  to  design  and  man¬ 
ufacture  more  efficient  and  powerful  gas  turbines  is  to  increase  the  turbine 
inlet  temperature,  for  the  power  output  and  thermal  efficiency  of  gas  turbine 
increases  with  the  increasing  inlet  temperature.  For  this  reason,  the  temper¬ 
ature  of  combustion  gas  at  turbine  inlet  at  present  far  exceeds  the  melting 
point  of  turbine  blade  superalloy,  and  very  sophisticated  measures  (e.g.  film 
cooling,  internal  cooling)  must  be  applied  to  cool  the  turbine  blades.  To  de¬ 
sign  effective  cooling  systems,  accurate  heat  transfer  prediction  between  the 
gas  and  the  turbine  blades  is  crucial.  Overestimating  heat  transfer  would 
cost  excessive  amount  of  cooling  air  bled  from  the  compressor,  hence  incur 
severe  penalties  in  the  engine  performance  and  efficiency.  Underestimating 
heat  transfer,  on  the  other  hand,  could  result  in  component  failure.  It  has 
been  reported  that  a  ten  degree  underestimation  of  turbine  blade  tempera¬ 
ture  would  reduce  the  blade  life  time  by  half.  Accurate  prediction  of  turbine 
heat  transfer,  however,  has  proven  to  be  rather  difficult.  One  of  the  key 
difficulties  stems  from  the  large  scale  and  high  intensity  turbulence  gener¬ 
ated  in  the  combustor  through  its  intricate  swirling,  mixing  and  combustion 
processes.  The  turbulence  contained  in  the  free  stream  greatly  enhances  the 
heat  transfer  between  the  combustion  gas  and  turbine  blades.  This  augmen¬ 
tation  effect  becomes  most  salient  in  the  stagnation  region  of  blade  leading 
edge,  where  the  largest  heat  transfer  occurs.  Because  the  physical  mecha¬ 
nism  which  leads  to  this  significant  heat  transfer  enhancement  is  not  well 
understood,  the  prediction  of  heat  transfer  in  the  presence  of  free-stream 
turbulence  remains  empirical  in  gas  turbine  design. 

Similar  heat  transfer  problems  related  to  stagnation  point  flow  also  arise 
in  micr<>electronics  cooling.  As  the  number  of  transistors  on  a  standard¬ 
sized  microprocessor  continues  to  increase,  the  power  density  is  expected 
to  exceed  100  W/cm^  in  the  near  future.  The  thermal  management  of  the 
power  dissipation  has  become  one  of  the  major  limiting  factors  in  building 
more  powerful,  compact  and  reliable  computers.  With  the  development  of 
MEMS  technology,  microjet  (or  jet  array)  impingement,  capable  of  producing 
heat  transfer  coefficient  one  order  of  magnitude  higher  than  that  of  the  tra¬ 
ditional  fan-driven  convection,  has  emerged  as  a  promising  new  technology. 
Active  experimental  researches  are  conducted  to  explore  various  methods 
for  increasing  the  heat  removal  rate  as  well  as  enhancing  the  controllability. 
Obviously,  a  detailed  knowledge  of  the  impinging  jet  flow  and  surface  heat 
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transfer  is  essential  to  achieving  a  optimized  cooling  design  which  ensures 
both  the  performance  and  reliability  for  future  microprocessors. 

Many  other  engineering  applications,  such  as  combustion  spray,  fluidized 
bed,  paper  drying,  wind-building  interaction,  need  to  deal  with  problems  in¬ 
volving  stagnation  point  flow  and  heat  transfer.  In  all  these  practical  cases, 
disturbances  exist  in  the  free  stream  and  in  fact,  are  often  the  major  concerns 
of  the  problem.  Hence,  a  better  understanding  of  the  dynamics  of  the  stag¬ 
nation  flow  and  heat  transfer  under  various  disturbances  is  at  the  heart  of  the 
fluid  mechanics  research  which  aims  at  improving  the  existing  technologies, 
and  is  also  the  primary  goal  of  the  present  study. 

In  what  follows,  past  studies  through  experimental  measurements,  the¬ 
oretical  analysis  and  numerical  computations  on  the  stagnation  point  flow 
and  heat  transfer  are  reviewed. 


1.2  Experimental  Studies 

The  experimental  studies  of  stagnation  point  flow  in  the  presence  of  distur¬ 
bances  started  in  later  1920s  when  the  unexpectedly  high  velocity  fluctuations 
were  first  observed  in  the  stagnation  region  of  a  blunt  body  placed  in  the  wind 
tunnels  with  uniform  streams  (Piercy  and  PUchardson,  1928,  1930).  Later, 
large  disparities  in  the  related  heat  transfer  measurements  among  different 
researchers  led  to  the  discovery  that  much  of  the  difference  can  be  attributed 
to  the  different  levels  of  the  free-stream  turbulence  in  the  wind  tunnels.  Fur¬ 
ther  systematic  experimental  studies  confirmed  that  the  free-stream  turbu¬ 
lence  causes  significant  heat  transfer  enhancement  in  the  stagnation  point 
region  (Giedt,  1949;  Hegge-Zijnen,  1957;  Kestin  et  al.,  1961). 

To  understand  and  model  the  effect  of  the  free-stream  turbulence  in  stag¬ 
nation  point  flows,  Smith  and  Kuethe  (1966)  used  eddy  viscosity  which  was 
assumed  to  be  proportional  to  the  product  of  free-stream  turbulence  intensity 
Tu  and  the  distance  from  the  wall.  Upon  applying  the  Hiemenz  transforma¬ 
tion,  the  combination  of  emerges  as  a  parameter  that  controls  the  shape 
of  the  turbulent  mean  velocity  and  thermal  profiles.  This  parameter  later 
forms  the  basis  of  many  empirical  correlations,  including  Smith  and  Kuethe 
(1966);  Kestin  and  Wood  (1971);  Lowery  and  Vachon  (1975);  Mehendale  et  al. 
(1991).  These  correlations  fit  the  FVossling  number  Fr  =  (Prossling, 

1940)  all  using  a  second  order  polynomial  of  but  with  slightly  different 
coeflScients  by  different  investigators. 
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Correlations  based  on  the  parameter  ^  achieved  some  success  in  data 
reduction  and  some  of  them  are  still  used  today.  But  they  generally  have  a 
rather  narrow  applicable  range  in  terms  of  Reynolds  number  and  turbulence 
intensity,  which  indicates  that  a  potentially  important  role  of  turbulence  scale 
has  been  missing  in  the  formulation.  The  reason  that  in  early  experiments 
the  length  scale  of  turbulence  did  not  receive  much  attention  —  most  of  the 
investigators  did  not  measure  the  length  scale  except  Hegge-Zijnen  (1957), 
are  probably  twofold.  In  the  first  place,  it  is  not  very  easy  to  achieve  experi¬ 
mentally  a  large  enough  turbulence  scale  range  so  that  a  definite  conclusion 
of  its  effect  can  be  drawn.  Secondly,  it  was  expected  that  the  turbulence 
length  scale  only  plays  a  subordinate  role  (Kestin  and  Wood,  1971).  As  was 
then  understood,  the  small  velocity  fluctuations  in  the  free  stream  triggers 
the  instability  of  the  mean  flow,  but  the  resulting  structures,  presumably  re¬ 
sponsible  for  the  heat  transfer  enhancement,  are  determined  by  the  intrinsic 
instability  features  rather  than  by  the  free-stream  turbulence  per  se. 

Presently,  the  importance  of  turbulence  length  scales  has  been  well  rec¬ 
ognized.  Intuitively,  very  small  scale  disturbances  are  quickly  damped  by 
the  viscosity,  whereas  very  large  scale  disturbances  only  exert  a  quasi-steady 
effect.  Hence  there  exists  a  certain  length  scales  range  in  which  the  distur¬ 
bances  are  most  effective.  However,  it  is  not  clear  what  is  the  precise  role 
of  the  turbulence  length  scales  in  determining  the  overall  turbulence  effect, 
and  different  opinions  exist  on  even  which  turbulence  scale  is  relevant  and 
should  be  incorporated  into  the  formulation. 

Yardi  and  Sukhatme  (1978)  measured  the  effect  of  the  turbulence  integral 
length  scale  L  on  the  stagnation  point  heat  transfer  using  a  circular  cylinder 
in  a  cross-flow,  and  showed  that,  at  both  very  large  and  very  small  scales, 
the  heat  transfer  enhancement  is  reduced.  The  maximum  effect  of  turbulence 
length  scale  appeared  somewhere  between  5  <  <15,  where  D  is  the 

diameter  of  the  cylinder,  or  roughly  4-12  times  of  boundary  layer  thickness. 

Ames  and  Moffat  (1990)  studied  the  effect  of  free-stream  turbulence  with 
relative  large  scales  (^  >  1.0)  on  stagnation  point  heat  transfer.  The  eddy 
viscosity  formulation  in  Smith  and  Kuethe  (1966)  was  generalized  by  tak¬ 
ing  into  account  the  change  of  the  turbulence  intensity  as  the  stagnation 
point  is  approached.  Based  on  a  model  turbulence  energy  spectrum  and 
rapid  distortion  theory  (Hunt,  1973;  Hunt  and  Graham,  1978),  a  new  correla¬ 
tion  parameter  is  developed  incorporating  the  turbulence  intensity,  Reynolds 
number  and  turbulence  energy  scale  (the  dissipation  scale  of  Hancock  and 
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Bradshaw  (1983))  in  the  form  of  The  energy  scale  represents  the 

energy-containing  eddies  and  is  estimated  by  Lu  =  where  c  is  the 

turbulence  energy  dissipation  rate.  This  parameter  appears  to  well  corre¬ 
late  their  experimental  data.  Data  from  other  groups  are  also  correlated  but 
with  more  scatter.  Further  experiments  have  been  reported  recently  (Ames 
et  al.,  2002)  with  a  larger  Reynolds  number  range  and  different  turbulence 
conditions. 

Van  Fossen  et  al.  (1995)  made  a  systematic  study  on  the  influence  of 
various  flow  parameters  on  stagnation  point  heat  transfer,  focusing  on  the 
relatively  short  turbulence  length  scale  (-^  =  0.05  —  0.3).  The  turbulence 
intensity  Tu  varies  from  1.1  to  15.9  percent,  and  Reynolds  number  Rep  from 
37000  to  228000.  They  showed  that  the  heat  transfer  enhancement  increases 
as  turbulence  length  scale  decreases,  but  no  peak  was  found  within  their 
range  of  length  scale.  The  measured  heat  transfer  enhancements  are  found 
to  be  best  fitted  using  an  empirical  parameter  of  Different  mean 

flow  strain  rates  are  also  studied  but  no  significant  influence  was  observed  on 
the  heat  transfer  enhancement.  Further  experiments  concerning  the  effect  of 
anisotropy  of  the  turbulence  (Oo  and  Ching,  2001),  and  for  even  higher  (upto 
28.5  percent)  turbulence  intensities  (VanFossen  and  Bunker,  2000)  have  also 
been  reported. 

Dullenkopf  and  Mayle  (1995)  argued  that  across  the  entire  turbulence 
energy  spectrum,  only  a  rather  selective  frequency  band  contributes  signif¬ 
icantly  to  the  heat  transfer  enhancement.  Thus,  instead  of  using  the  total 
turbulence  intensity,  they  proposed  to  use  an  effective  turbulence  intensity 
obtained  by  integrating  the  von  Karman’s  model  spectrum  over  a  narrow 
frequency  range.  Based  on  their  own  and  some  previous  data,  they  found 
that  the  heat  transfer  enhancement  appears  to  vary  linearly  with  this  length 
scale  weighted  turbulence  intensity,  reminiscent  to  the  earlier  prediction  by 
Smith  and  Kuethe  (1966). 

It  should  be  pointed  out  that,  in  most  cases,  correlations  developed 
through  specific  experiments  are  not  applicable  to  data  from  other  researchers. 
Moreover,  of  particular  interest  to  the  gas  turbine  applications  is  the  stag¬ 
nation  point  flow  and  heat  transfer  with  turbulence  present  in  a  transonic 
free-stream,  but  experimentally  this  has  proven  very  challenging  and  few 
measurements  are  available.  Difficulties  also  arise  when  the  correlations  de¬ 
veloped  using  simplified  geometry  are  to  fit  the  heat  transfer  data  in  a  turbine 
blade  leading  edge  region  (Yeh  et  al.,  1993;  Thole  et  al.,  2002).  On  one  hand. 
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this  may  be  attributed  to  the  sensitivity  of  the  stagnation  point  flow  to  var¬ 
ious  free  stream  disturbances.  On  the  other  hand,  the  empirical  nature  of 
these  experimental  correlations  is  the  inherent  reason  for  their  limited  appli¬ 
cability.  From  this  perspective,  theoretical  analysis  should  shed  more  light 
on  the  underlying  physical  mechanism,  and  lead  to  a  better  understanding 
of  the  precise  roles  of  the  various  flow  parameters  used  in  the  experimental 
correlations. 


1.3  Theoretical  Analysis 

In  an  attempt  to  explain  the  sensitivity  of  the  stagnation  point  flow  and  heat 
transfer  to  the  upstream  disturbances,  theoretical  analysis  has  been  carried 
out  from  various  perspectives,  mainly:  unsteady  mean  flow  effect,  linear  and 
nonlinear  instability  and  vortex  stretching  and  ampliflcation. 

The  unsteadiness  of  the  mean  flow  caused  by  upstream  disturbances  was 
first  considered  by  Lighthill  (1954),  who  examined  the  effect  of  pulsating  up¬ 
stream  mean  velocity  about  a  cylindrical  body  and  obtained  the  Stokes-layer 
correction,  but  no  significant  change  in  heat  transfer  was  found.  Similarly 
modulated  stagnation  point  flows  were  also  studied  by  Ishigaki  (1970)  and 
Pedley  (1972),  but  with  an  emphasis  on  the  surface  skin  friction.  Merchant 
and  Davis  (1989)  examined  the  steady  streaming  driven  by  large  amplitude 
and  high  frequency  temporal  free  stream  modulation.  The  results  show  that 
the  streaming,  when  it  exists,  is  always  confined  within  the  Hiemenz  bound¬ 
ary  layer  thickness.  These  unsteady  analyses  are  limited  to  two  dimensional 
case,  thus  do  not  explain  the  large  three  dimensional  structures  often  ob¬ 
served  in  the  perturbed  stagnation  point  flow.  These  structures,  presumably 
the  reason  of  the  heat  transfer  enhancement,  often  take  the  form  of  an  array 
of  secondary  streamwise  vortices  separated  by  a  distances  of  the  order  of 
boundary  layer  thickness. 

From  the  point  of  view  of  flow  instability,  these  structures  are  hypoth¬ 
esized  as  the  result  of  a  selective  amplification  process,  i.e.  the  three  di¬ 
mensional  disturbances  originating  from  the  viscous  stagnation  region  of 
one  particular  spanwise  wavelength  is  preferentially  amplified.  The  wave¬ 
length  may  potentially  correspond  to  the  distance  of  vortex  spacing  observed 
in  the  experiments.  Early  linear  stability  analysis  by  Gortler  (1955)  and 
Hammerlin  (1955),  however,  showed  that  the  resulting  eigenvalue  problem 
yields  a  continuous  spectrum  of  spanwise  wavenumbers  for  the  neutral  time 
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dependent  disturbance  as  opposed  to  a  unique  most  amplified  wave  number. 
Later,  Kestin  and  Wood  (1970)  pointed  out  that  the  lack  of  unique  eigen¬ 
value  is  the  result  of  the  over  simplification  of  the  boundary  conditions  at 
far  upstream.  The  mathematical  formulation  was  considerably  clarified  by 
Wilson  and  Gladwell  (1978),  who  showed  that  the  correct  behavior  of  the 
disturbances  originating  from  the  stagnation  region  should  be  exponentially 
decaying  far  upstream.  The  resulting  stability  analysis  in  this  case  led  to 
an  unique  eigenvalue  problem,  but  the  disturbances  are  all  found  to  decay, 
i.e.  the  plane  stagnation  point  fiow  is  always  linearly  stable  to  the  three 
dimensional  disturbances.  The  nonlinear  instability  of  finite  amplitude  dis¬ 
turbances  was  further  considered  by  Lyell  and  Huerre  (1985)  using  Galerkin 
method,  and  they  showed  that  the  flow  can  be  destabilized  if  the  level  of  the 
external  two-  or  three  dimensional  disturbances  exceeds  certain  threshold 
value. 

In  a  comprehensive  review,  Morkovin  (1979)  argued  that  the  enhancement 
of  heat  transfer  is  more  likely  a  result  of  forced  response  to  the  upstream  dis¬ 
turbances  as  opposed  to  the  internal  flow  instability.  The  flow  visualizations 
by  Nagib  and  Hodson  (1978)  and  Botcher  and  Wedemeyer  (1989)  strongly 
support  this  argument.  The  same  point  of  view  was  advocated  earlier  by 
Sutera  (1965)  who  analyzed  the  amplification  effect  of  the  mean  flow  strain¬ 
ing  on  the  incoming  organized  disturbances,  indicating  the  sensitivity  of  the 
heat  transfer  to  vortical  disturbances.  Hunt  (1973)  used  rapid  distortion 
theory  to  investigate  the  asymptotic  behavior  of  free-stream  turbulence  as  it 
approaches  a  bluff  body;  considerable  insight  was  gained  into  the  distortion 
process  of  the  turbulence  of  different  length  scales.  But  the  mean  flow  was 
assumed  potential,  therefore  no  heat  transfer  between  the  fluid  and  wall  was 
considered.  To  explain  the  formation  of  the  secondary  vortices,  Kerr  and 
Dold  (1994)  showed  that  an  inviscid  two  dimensional  stagnation  flow  can 
posses  periodic  streamwise  vortices.  The  characteristics  of  these  vortices  are 
also  studied  in  the  experiments  of  Andreotti  et  al.  (2001).  To  connect  the  dis¬ 
turbances  in  the  free  stream  with  the  vortices  found  in  the  stagnation  region, 
Dhanak  and  Stuart  (1995)  showed  that  the  weak  external  cross-stream  vortic- 
ity  with  small  length  scales  follows  an  algebraic  structure  in  the  wall  normal 
direction,  under  which  the  corresponding  inner  viscous  boundary  layer  can 
support  a  substructure  of  counter-rotating  streamwise  eddies.  Regarding  the 
effect  of  the  three  dimensional  flow  structures  on  the  heat  transfer,  the  nu¬ 
merical  simulation  of  Rigby  and  VanFossen  (1992)  showed  that  the  small 
spanwise  nonuniformity  can  result  in  large  increase  of  heat  transfer  in  the 
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stagnation  point  region  of  a  cylinder.  Bae  et  al.  (2000)  also  found  that  differ¬ 
ent  length  scales  generate  quite  different  flow  patterns  and  in  turn  different 
heat  transfer  responses  in  a  plane  stagnation  point  flow. 

Despite  the  significant  progress  made  in  the  past,  the  precise  mechanism 
governing  the  disturbance  evolution  and  heat  transfer  augmentation  is  still 
not  fully  understood,  and  particularly  elusive  is  its  dependence  on  the  various 
disturbance  parameters. 


1.4  Numerical  Predication 

To  represent  the  turbulence  effect  on  the  flow  and  heat  transfer  for  most  of  the 
practical  applications,  the  Reynolds  averaged  Navier-Stokes  equations  based 
on  various  turbulence  models  are  at  present  still  the  major  predication  tools. 
However,  stagnating  turbulence  presents  a  serious  challenge  for  turbulence 
modelling.  This  is  chiefly  because  the  formulation  of  turbulence  models  is  in 
most  cases  based  on  the  shear  dominated  flows  such  as  boundary  layer  rather 
than  the  strain  dominated  flows  like  those  found  in  a  stagnation  point. 

One  equation  model  of  Spalart  and  Allmaras  (1992)  has  been  widely 
used  and  shown  particularly  successful  in  aerodynamic  flows  (Bardina  et  al., 
1997;  Wilcox,  2001).  In  this  model  the  eddy  viscosity  is  obtained  algebraically 
through  an  effective  eddy  viscosity  which  satisfies  a  proposed  transport  equa¬ 
tion.  The  model,  however,  does  not  explicitly  account  for  the  effect  of  free- 
stream  turbulence.  In  order  to  predict  the  effect  of  free-stream  turbulence, 
e.g.  heat  transfer  enhancement,  bypass  transition,  modifications  must  be  in¬ 
troduced  into  the  formulation  (Tsourakis  et  al.,  2002).  Those  modifications 
are  often  of  rather  ad  hoc  nature,  and  it  is  not  clear  how  in  general  that  the 
effect  of  free-stream  turbulence  may  be  incorporated. 

Standard  two  equation  models,  e.g.  k  -  e  or  k  -  u  model,  when  used 
in  the  stagnation  point  turbulent  flows,  badly  over  predict  the  turbulent  ki¬ 
netic  energy  and  heat  transfer  —  “stagnation  anomaly”  has  been  termed  to 
describe  this  phenomenon  (Durbin,  1996)  although  it  happens  wherever  the 
mean  strain  rate  is  large  in  the  flow  field,  such  as  in  a  gas  turbine  passage 
(Medic  and  Durbin,  2002).  For  uniform  mean  flow  strain,  rapid  distortion 
theory  shows  that  the  production  of  the  turbulent  kinetic  energy  is  pro¬ 
portional  to  the  mean  flow  strain  rate.  But  in  the  two  equation  models, 
the  normal  stresses  are  expressed  by  eddy  viscosity  and  mean  flow  gradi¬ 
ents,  thus  the  production  terms  becomes  proportional  to  the  square  of  the 
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strain  rate.  Moreover,  turbulence  effect  in  two  equation  models  is  repre¬ 
sented  chiefly  by  its  kinetic  energy,  implying  the  turbulence  is  isotropic.  But 
in  stagnation  point  flows,  the  free-stream  turbulence  quickly  becomes  very 
anisotropic,  because  different  components  of  the  velocity  fluctuation  respond 
differently  to  the  mean  flow  strain  rate  when  the  wall  is  approached.  These 
two  factors  essentially  render  the  standard  two  equation  model  inapplicable 
in  the  stagnation  point  type  of  flows.  Various  remedies  have  been  proposed 
to  cure  the  problem,  ranging  from  modifying  the  standard  model  coefficients 
(Champion  and  Libby,  1991),  changing  the  upstream  boundary  condition 
(Abid  and  Speziale,  1994),  reformulating  the  turbulence  production  term 
(Kato  and  Launder,  1993)  to  zonal  matching  (Traci  and  Wilcox,  1975;  Wang 
and  Yeh,  1987)  and  using  second  order  Reynolds  stress  model  for  the  nor¬ 
mal  stresses  (Taulbee  and  Tran,  1988).  One  recent  development  is  to  use 
time-scale  limiters  (Durbin,  1996),  where  an  additional  bound  for  the  turbu¬ 
lent  characteristic  time  scale  is  derived  from  the  realizability  considerations. 
When  implemented  into  standard  two  equation  models,  it  shows  promising 
results  in  predicting  the  correct  turbulence  behavior  and  heat  transfer  coef¬ 
ficient  (Medic  and  Durbin,  2002;  Behnia  et  al.,  1999;  Parneix  et  al.,  1999). 

The  Reynolds  stress  models  are  the  most  complicated  turbulence  models. 
They  are  used  less  often  than  the  two  equation  models  because  not  only 
more  equations  are  involved,  but  also  the  coupling  among  different  Reynolds 
stress  terms  makes  the  numerical  problem  difficult.  Nevertheless,  in  Reynolds 
stress  models,  the  anisotropy  of  the  turbulence  is  explicitly  accounted,  so 
better  predicative  capability  may  be  expected.  However,  though  better  than 
k  —  €  type  models,  Reynolds  stress  models  are  still  far  from  satisfactory  when 
applied  to  stagnating  turbulent  flows  such  as  an  impinging  jet.  Im  et  al. 
(2003)  used  three  variants  of  Reynolds  stress  model,  GL  model  (Gibson  and 
Launder,  1978),  GL-CL  model  (Craft  et  al.,  1993)  and  SSG  model  (Speziale 
et  al.,  1991),  to  compute  both  the  impinging  and  counter-current  stagnation 
flows.  All  the  models  give  severe  over-predictions  in  turbulence  kinetic  energy 
and  large  discrepancies  in  other  Reynolds  stress  components  when  compared 
with  experimental  measurements.  The  problem  stems  from,  on  one  hand  the 
over-prediction  of  the  energy  production,  while  on  the  other  hand  the  under 
prediction  of  the  redistribution  due  to  the  troublesome  terms  modelling  the 
pressure  strain  correlations. 

In  summary,  from  the  aforementioned  studies  using  experimental,  the¬ 
oretical  and  numerical  approaches,  it  is  clear  that  despite  the  effort  and 
progress  made  over  the  years,  many  important  questions  regarding  to  the 
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perturbed  stagnation  point  flow  remain  to  be  answered.  A  closer  and  more 
detailed  investigation,  combining  both  numerical  simulation  and  theoretical 
analysis,  is  thus  of  great  interest.  Not  only  will  it  shed  much  light  into  the 
physical  mechanism  but  also  can  potentially  provide  valuable  guidance  to¬ 
wards  developing  and  calibrating  new  turbulence  models.  The  insight  gained 
from  the  study  may  improve  the  design  of  experiments  and  help  developing 
new  technologies  in  a  wide  range  of  related  industrial  fields. 


1.5  Overview  and  accomplishments 

In  this  study,  we  adopt  both  numerical  and  theoretical  approaches  to  inves¬ 
tigate  the  the  turbulent  flow  and  heat  transfer  at  a  two  dimensional  elliptical 
leading  edge.  The  principal  contributions  and  findings  of  this  work  are  listed 
below. 

•  High  order  fully  implicit  numerical  method,  with  a  linearized  dual  time 
sub-iteration  scheme  offering  a  4-5  times  speed-up  over  the  usual  subit¬ 
eration  schemes,  is  developed  and  validated  to  study  the  stagnation 
point  flow  and  heat  transfer  in  the  presence  of  free-stream  turbulence. 

•  A  new  blending  procedure  is  developed  for  generating  realistic  free 
stream  turbulence.  By  blending  a  series  of  independent  but  statistically 
identical  realizations  of  homogeneous  isotropic  turbulence,  this  simple 
yet  effective  method  preserves  the  first  order  and  most  of  the  second 
order  turbulence  statistics. 

•  For  low  Mach  number  case,  free-stream  turbulence  conditions  with  dif¬ 
ferent  length  scales  and  turbulence  intensities  matching  experimental 
conditions  are  computed  using  large  eddy  simulation.  Very  good  agree¬ 
ment  between  the  numerical  simulation  and  experimental  measurement 
on  the  heat  transfer  enhancement  (Van  Fossen  et  al.,  1995)  is  obtained. 

•  Three  different  stages  characterizing  the  interaction  between  free-stream 
turbulence  and  stagnation  point  are  identified  and  examined  in  detail. 
Turbulence  statistics  are  obtained  and  analyzed  in  the  light  of  rapid 
distortion  theory. 

•  The  impinging  and  lateral  movement  of  stretched  free  stream  turbu¬ 
lent  eddies  are  found  to  be  directly  responsible  for  the  heat  transfer 
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enhancement.  The  wall  blockage  efTect  are  shown  to  be  the  cause  of 
rapid  lateral  movement  of  the  eddies,  and  the  magnitude  of  the  trans¬ 
lating  velocity  is  characterized. 

•  Large  eddy  simulation  at  high  Mach  number  is  performed,  and  the 
turbulence  statistics  are  obtained  and  compared  with  the  low  Mach 
number  case. 

•  The  length  scale  effect  on  the  disturbance  evolution  is  analyzed  us¬ 
ing  matched  asymptotic  expansion.  A  new  uniformly  valid  analytic 
expression  is  obtained  which  describes  evolution  of  three  dimensional 
unsteady  upstream  disturbances  being  convected  towards  the  stagna¬ 
tion  point.  The  effect  of  mean  flow  straining  and  viscous  dissipation  as 
well  as  its  implication  on  heat  transfer  are  explicitly  shown. 

•  The  analysis  for  organized  disturbances  is  generalized  to  the  free  stream 
turbulence.  A  new  correlation  parameter  is  developed  which  corre¬ 
lates  the  heat  transfer  enhancement  to  the  turbulence  intensity,  inte¬ 
gral  length  scale  and  Ileynolds  number.  The  new  parameter  is  shown  to 
reasonably  collapse  both  the  recent  experimental  data  and  the  results 
of  the  present  numerical  simulations. 

This  report  documents  the  primary  computation  and  analysis  methods  as 
well  as  the  results  and  findings  at  the  time  of  writing.  More  detailed  technical 
information  will  be  available  in  the  Stanford  University  Ph.D  thesis  (Xiong, 
2004). 
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Chapter  2 

Governing  Equations 


2.1  Governing  Equations  for  Compressible  Flow 

The  physical  laws  governing  the  motion  of  a  compressible  fluid  are  expressed 
through  the  continuity,  momentum  and  energy  equations.  Using  Cartesian 
tensor  notation  and  dimensional  primitive  variables,  they  can  be  written  as 

Pit  +  (/>*<),<  =  0  (2.1) 

=  -Pli  +  Ttjj  (2.2) 

p*h*t  +  p*u*jh*j  =  {p*t  +  «>*-)  +  qjj  +  (2.3) 

where  p*  is  the  density,  is  the  velocity  vector,  p*  is  the  thermodynamic 
pressure,  is  the  viscous  stress  tensor,  Qj  is  the  heat  flux,  and  h*  is  the  fluid 
enthalpy  defined  by 

h*  =  e*+p*/p*  (2.4) 

where  e*  is  the  internal  energy  per  unit  mass.  In  the  above  equations,  sub¬ 
scripts  following  a  comma  denote  partial  differentiation  with  respect  to  the 
subscript,  and  the  Einstein  summation  convention  is  used. 

We  assume  Newtonian  fluid  such  that  the  viscous  stress  tensor  and  the 
heat  flux  are  given  by 

(2.6) 

«;  =  --'Tj  (2.6) 

where  p*  and  A*  are  the  first  and  second  coefficient  of  viscosity,  k*  is  the 
thermal  conductivity,  =  2(^ij  +  ^i,t)  rate-of-strain  tensor,  T*  is  the 

absolute  temperature,  and  %  is  the  Kronecker  delta. 
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To  close  the  system  of  equations  (2.1)  -  (2.3),  an  equation  of  state  which 
links  the  thermodynamical  variables  and  additional  equations  which  relate 
the  fluid  properties  to  the  thermodynamical  variables  are  needed.  We  assume 
an  ideal,  calorically  perfect  gas  with  constant  specific  heats,  so  the  equation 
of  state  is 

p*  =  p*n*T*  (2.7) 

were  H*  =  c*  —  c*  is  the  gas  constant,  c*  is  the  specific  heat  at  constant 
pressure,  and  c*  is  the  specific  heat  at  constant  volume.  Both  c*  and  c* 
are  constant  as  is  the  ratio  7  =  c*/c*.  Under  these  conditions,  the^internal 
energy  and  the  enthalpy  are  related  to  the  absolute  temperature  by 

e*  =  <T*  (2.8) 

=  c;T*  (2.9) 

The  fluid  properties,  for  a  calorically  perfect  gas,  are  functions  of  temperature 
only.  For  low  speed  flows  with  small  temperature  variations,  the  viscosity 
p*  and  conductivity  /c*  can  be  treated  as  constants.  When  the  temperature 
dependence  becomes  important,  p*  and  k*  are  often  expressed  in  empirical 
relations  in  a  nondimensional  form  involving  reference  states.  So  before  de¬ 
scribing  them,  we  first  nondimensionalize  (2.1)  -  (2.9)  using  the  following 
expressions. 


p  = 

£l 

Pr' 

“«■  = 

rp  _  T* 

■“  r;  ’ 

T%  - 

V 

P 

A  —  — r, 

K  = 

K* 

~  ’ 

t  = 

where  subscript  r  denotes  the  reference  variables,  whose  particular  values 
are  defined  for  each  problem  considered.  Using  (2.10)  in  (2.1)  -  (2.7),  the 
nondimensional  governing  equations  take  the  form 

P,t  +  (pMi),i  =  0  (2.11) 

pUi,t  +  pUjUij  =  -pi  -f-  —[[Xuj  .)  .  +  {2pSij)j]  (2.12) 

pTt-^pUjTj  -I-  (7-l)pTujj  =  -^[KTi],i-i- 

7(7- 1)^2,, ^  ^ 

Re  (2.13) 
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where 


„  PT 
^  7M2 


(2.14) 


Re  = 


P>*rL; 

Pr 


Pr=^ 


(2.15) 


are  the  Mach  number,  Reynolds  number  and  Prandtl  number.  In  deriving  the 
nondimensional  energy  equation  (2.13),  (2.9)  is  used  to  convert  enthalpy  to 
temperature.  The  continuity  equation  (2.1)  along  with  the  equation  of  state 
(2.7)  are  used  to  remove  the  explicit  pressure  dependence.  The  nondimen¬ 
sional  viscosity  //  =  1  when  it  is  independent  of  temperature,  and  otherwise 
a  power  law  is  used 

/i  =  T”  (2.16) 


where  n  is  taken  as  n  =  0.70.  The  second  coeflScient  of  viscosity  is  computed 
using  Stokes  hypothesis,  A  =  — |/i  (giving  zero  bulk  viscosity),  and  the  ther¬ 
mal  conductivity  is  determined  by  assuming  constant  Prandtl  number,  so 
that  K  =  fx.  Unless  otherwise  specified,  the  fluid  is  assumed  to  be  air  with  a 
Prandtl  number  of  0.71  and  a  ratio  of  specific  heats  j  —  1.4. 


2.2  Filtered  Governing  Equations 

The  governing  equations  for  large  eddy  simulation  (LES)  are  obtained  by 
filtering  equations  (2.11)-(2.14).  Let  a  filtered  or  large-scale  flow  quantity  be 
denoted  by  an  overbar 


f{x)=  f  G{x  —  3/)f{x')dx'  (2.17) 

Jd 

where  G  is  some  spatial  filter  and  the  integral  is  over  the  flow  domain  D. 
The  velocity  tii,  density  p  and  temperature  T  can  be  decomposed  as 

Ui  =  Ui p  =  p (J ,  T  =  T -\-T'  (2.18) 

Applying  the  spatial  filter  G  to  equation  (2.11)-(2.14)  leads  to 

P,t  +  =  0  (2.19) 

=  -^  +  +  {2pSij)J\  (2.20) 
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(2.21) 


pTt  +  pUiTi  +  {^-\)pTuij 

7(7  -  1)M2 


PrRe 


4-  2pSijSij 


P  = 


pT 

7M2 


(2.22) 


The  large-scale  field  equations  (2.19)-(2.22)  can  be  operationally  simplified 
when  the  variables  are  recast  in  terms  of  Favre-filtered  quantities.  A  Farve- 
filtered  variable  is  defined  as 


f  =  pf/p  (2.23) 

Note  that 

Wi  =  pUi,  puj^  =  pi^,j,  ^jTj  =  p^j  (2.24) 

The  stress  tensor  and  the  heat  flux  can  be  decomposed  into  resolved  and 
sub-grid  scale  (SGS)  components  as  follows: 

pupTj  =  'pUiUj  4-  'p{u^j  -  UiUj)  (2.25) 

> - - -  / 

'pUiT  =  puif  +  p{uiT  -  Ujf)  (2.26) 

Similar  decompositions  are  made  for  the  viscous  terms  in  the  momentum  and 
energy  equations  and  the  pressure  dilatation  and  conduction  terms  in  energy 
equation.  But  their  small  scale  subgrid  components  were  neglected.  When 
the  filtering  operation  (2.17)  is  applied  to  homogenous  directions  of  the  flow, 
it  commutes  with  the  differentiation  operator.  In  this  case,  the  governing 
equations  for  the  large-scale  field  become 

P,t  +  ^i),i  =  0  (2.27) 

pui,i  +  pUjUij  =  -Pi  4-  +  (2/i4).,]  -  T^. (2.28) 

pf.  +  piif,  +  {l-DWuij=p^[kfA,+ 

7(7  -  l)A/2  -  .  -  _  . 

Re  ,  (2,29) 
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(2.30) 


^  7M2 

where  fi  is  the  same  as  in  (2.16)  except  the  T  is  replaced  by  T.  The  second 
viscous  coefficient  A  is  still  given  by  A  =  -|/i  and  thermal  conductivity 

K  =  K. 


2.3  Dynamic  Model  for  SGS  Transport 

To  close  the  momentum  and  energy  equation,  and  must  be  modelled. 
For  notational  convenience,  the  '  will  be  omitted  in  what  follows.  For  r^,  the 
trace-free  Smagorinsky  eddy  viscosity  model  is  used 

nj  =  -  2C?A^|S|(4  -  (2.31) 


where  =  tu  is  the  isotropic  part  of  the  SGS  Reynolds  stress  tensor,  Sij  = 
and  |,S|  =  (25, The  SGS  energy  is  parametrized  using 
Yoshizawa’s  expression: 

q^  =  2C/pA2|5p  (2.32) 

For  the  qi,  the  eddy  diffusivity  model  is  used 


Qi  = 


C^A2|5| 


(2.33) 


where  the  Prt  is  the  SGS  turbulent  Prandtl  number  and  C  is  the  same  as 
in  (2.31).  We  will  use  dynamic  procedure  to  compute  the  eddy  coefficients 
C,  Cj  and  SGS  turbulent  Prandtl  number  Prt  (Moin  et  al.,  1991).  For  this 
purpose,  first  define  a  test  filter,  denoted  by  a  caret,  with  a  width  A  larger 
than  that  of  the  resolved  grid  filter.  The  test-filtered  stresses  J-ij  is  defined 
by  direct  analogy  to  the  stresses  : 

^ij  =  W^j  -  iW)  (^)  /  f  (2.34) 


By  Germano’s  identity  (Germano,  1990),  the  Leonard  stresses  Cij  can  be 
expressed  in  terms  of  and  as 

Cij  =  Tij  -  fij  =  (pUiUj)  -  (^f)  i'puj)  /  ^  (2.35) 
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Since  the  test  filter  is  always  applied  to  the  resolved  field,  all  the  quantities 
with  "are  obtained  from  the  computed  field.  So  the  right-hand  side  of  (2.35) 
is  computable  from  the  solved  variables.  To  determine  coefficient  Cj,  the 
trace  of  (2.35)  is  used  in  conjunction  with  the  model  of  (2.31)  for  Tu  and  th: 

'j^Ui  -  {pui)  (pUi)/^  =  2C'/(^A^|5p  -  A2p|5P)  (2.36) 

Hence  Cj  can  be  written  as 

C,  =  <  >  ^2  37) 

<  2(M^|5|2  -  A2pl5|2)  > 

where  <>  indicates  some  kind  of  volume  averaging  procedure  which  is  needed 
to  make  the  determination  of  Cj  and  the  other  SGS  coefficients  well  condi¬ 
tioned.  In  this  LES  study,  Cj  =  0  has  been  chosen  following  most  of  the 
previous  studies  for  numerical  stability. 

To  obtain  C,  we  also  use  Smagorinsky  model(2.31)  for  the  test  field 
stresses,  i.e. 


—  -^kkSii  = 


kk'^ij 


-2C^A2|||(4  -  ^SkkSij) 


(2.38) 


Substituting  and  into  (2.35),  after  contracting  with  Sij  and  appropri¬ 
ate  spatial  averaging,  we  obtain 

-  {'^i){'^j)/plSij  -  ISkkj^mm  -  fmm)  > _ 

<  -2^pA^\'s\CSijSij  -  fSkkSmm)  +  2A2(p'jJ|4-4.  -  0\§kkSmm)  > 

(2.39) 

where 

^ mm  ~  Tmm  —  P^m^m  ~  {p^m)  {p^m)/ P  (2.40) 

To  determine  Pn,  let  Qi  denote  the  the  heat  flux  at  the  test  filter  scale 
and  use  the  same  eddy  diffusivity  model,  it  follows  that 


(?,  =  -^t  -  (pUi)  (pf )  fp  =  (2.41) 

Substituting  Qi  and  qi  into(2.35),  contracting  with  Tj  and  performing  the 
appropriate  spatial  averaging,  we  obtain: 


C  <  A2^f,  -  A^|T,f,  > 

<  {pUi  ^  /p  -  ^if)fi  > 


(2.42) 
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where  C  is  given  by  (2.39). 

With  the  dynamic  SGS  model,  the  equations  for  the  filtered  quantities 
now  take  the  form 

P,t  +  =  0  (2-43) 

-  -PT,i  +  4-  (2/ir5ti)  j]  (2-44) 


7 


+  pUiT^i  +  (7  —  Pj-jiq 

7(7  -  1)M^ 


Re 


[AiSjjS'jj  +  2jiSijSij] 


(2.45) 


where 


^  |C/PA’|SP 

Pt  ~  P  +  (7pA^|iS|/?e 


kT-=k  +  CpA^\S\ 


X,  PrRe 


jPrt 


(2.46) 

(2.47) 

(2.48) 

(2.49) 


This  system  of  equations  constitutes  the  governing  equations  for  the  large 
eddy  simulation  of  compressible  turbulence  in  the  present  study. 
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Chapter  3 

Numerical  Method 


We  describe  in  this  chapter  the  numerical  method  of  solving  the  compress¬ 
ible  flow  governing  equations  (2.43)  -  (2.49).  The  equations  are  discretized 
using  a  fully  implicit,  approximately-factorized  (AF)  time  marching  scheme 
in  conjunction  with  fourth  order  spatial  center  differencing  in  a  curvilinear 
coordinate.  A  linearized  dual  time  sub-iteration  scheme,  which  offers  a  4  —  5 
times  speed-up  over  the  usual  subiteration  scheme,  is  developed.  The  details 
of  the  finite  difference  scheme,  implicit  boundary  conditions  and  other  as¬ 
pects  such  as  numerical  dissipation  are  also  discussed.  The  numerical  code  is 
validated  against  the  analytical,  self-similar  solution  of  stagnation  point  com¬ 
pressible  boundary  layer,  and  the  computation  of  leading  edge  receptivity  to 
sound  of  a  flat  plate  compressible  boundary  layer. 

3.1  Implicit  Scheme  with  Linearized  Subiter¬ 
ations 

Implicit  methods  have  been  traditionally  developed  to  compute  steady  state 
or  slowly  varying  unsteady  flows,  whereas  explicit  methods  like  R-K  schemes 
were  typically  the  choices  for  time  accurate  computations  in  the  past.  With 
the  growing  demand  for  detailed  flow  simulations  with  complicated  geometry 
particularly  in  the  presence  of  solid  surfaces,  the  high  grid  resolution  inside 
boundary  layers  often  makes  the  unsteady  implicit  method  more  preferable, 
because  larger  time  steps  can  be  taken  than  would  be  permitted  by  explicit 
stability  limits. 

For  the  implicit  time  marching  scheme,  we  first  recast  the  compressible 
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flow  governing  equations  into  a  general  form: 


j7,t  +  F(t/)  =  0  (3.1) 

where  U=  {p,  u,  v,  w,  T}^  is  the  vector  of  flow  variables  and  F{U)  represents 
the  nonlinear  and  viscous  terms.  With  the  second  order  backward  Euler 
scheme,  (3.1)  can  be  written  in  a  semi-discretized  form  as 


-  4f7«  +  ijn-i 

2At 


-f  F(f/"+i)  =  0 


(3.2) 


In  the  dual  time  formulation,  a  new  pseudo  time-derivative  is  intro¬ 
duced  into  (3.2),  i.e. 


Ur  + 


3f/n+l  -  4Un  ^  ijn-1 

2At 


-I-  F(C/"+i)  =  0 


(3.3) 


where  t  is  the  pseudo  time  variable.  Hence  when  (3.3)  reach  a  steady  state, 
i.e.  it  recovers  (3.2).  To  this  purpose,  we  first  linearize  F(f7«+i)  with 
respect  to  f/" 


F(Lr"+i)  =  F{ir)  +  A{ir)  Air  +  0{AU^)  (3.4) 

where  .4(17)  =  Fu  is  the  Jacobian  of  F{U)  and  AIT  =  Typically, 

in  a  multi-dimension  system,  the  size  of  the  coefficient  matrix  4,  resulting 
from  spatial  discretizations,  is  of  the  the  order  {NxNy  NzY,  and  often 
prohibitively  expensive  to  be  inverted  directly,  where  N^,  Ny,  are  number 
of  grid  points  in  x,  y,  z  directions.  In  standard  factorization  (AF)  method, 
the  matrix  4  is  replaced  by  a  sequence  easily  invertable  matrices  (Briley 
and  McDonald,  1975;  Beam  and  Warming,  1978)  whose  product  recovers 
approximately  to  4  with  factorization  errors  on  the  order  of  (A/)^.  To  ensure 
the  factorization  error  is  negligible  and  not  to  corrupt  the  solution,  a  time 
step  far  smaller  than  is  required  by  time  accuracy  has  to  be  used  in  unsteady 
computation. 

To  alleviate  this  problem,  an  inner  iteration  called  subiteration  is  intro¬ 
duced  for  each  physical  time  step  of  an  unsteady  computation  (Rai,  1987; 
Pulliam,  1993).  If  the  subiteration  converges,  the  factorization  error  will 
be  eliminated.  For  the  steady  state  computation,  subiteration  enhances  the 
stability  and  robustness  of  the  AF  method,  whereas  for  unsteady  computa¬ 
tion,  it  will  recover  the  desirable  time  accuracy.  But  the  efficiency  of  the  AF 
method  now  strongly  depends  on  the  convergence  rate  of  the  subiteration. 


22 


So  introduce  a  subiteration  index  k  and  approximate  U^t  by  first  order 
backward  Euler  scheme, 


U,r  = 


IJk+l  _  ijk 

At 


(3.5) 


(3.3)  becomes 


+  TTT-  -  4t/”  +  ir-^)  + 

2At 

^rA{ir)Mr  =  -ATF{tr)  (3.6) 


where  the  At  is  an  appropriately  chosen  pseudo  time  step  for  subiteration 
and  At/*  =  —  U^.  Furthermore,  replacing  n  +  1  by  +  1  and  subtract 

from  both  sides  of  the  3.6,  we  obtain 


(/  +  — )  Af/*  +  ArAilT)  (C/‘+^  -ir)  = 

-^(3t/''  -  +  C/"-i)  -  AtF{U^)  (3.7) 


Finally,  notice  that 


jjk+l  _ljn^  ^^k  _ 


(3.8) 


and  expand  F{U’‘)  at  t/”,  the  subiteration  scheme  for  AU^  can  be  expressed 
as 

(3.9) 


|/  +  +  AT.4(i7”)]AE/‘  =  -AtK* 


where 


^  317*  -  4U^  + 


77*  = 


2At 


+  F(f/*) 


(3.10) 


For  each  physical  time  step,  take  t/”  as  the  initial  value  for  t/*  where  =  0  to 
start  the  subiteration.  When  the  subiteration  converges,  we  have  At/*  -4  0, 
Ijk+i  _  fjfc  -pijg  Qf  ijk+1  jjg  taken  as  17"+^  and  =  0 

recovers  the  second  order  fully  implicit  scheme  in  (3.2). 

Notice  that  the  left  hand  side  operator  in  (3.9)  is  only  a  function  of  IT*, 
whereas  in  the  standard  subiteration  schemes,  it  is  a  function  of  17*.  So  (3.9) 
is  actually  linear  for  the  subiteration  variable  A17*  (Venkateswaran  et  al., 
1997);  we  can  first  perform  the  LU  decomposition  of  the  coefficient  matrices 
and  store  the  factored  matrices  in  the  first  step  and  use  them  throughout 
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until  the  subiteration  converges.  This  obviates  the  need  of  computing  Jaco- 
bians  and  inverting  the  coefficient  matrices  at  each  subiteration  step,  there¬ 
fore  significantly  improves  the  subiteration  efficiency.  Moreover,  the  present 
subiteration  formulation  is  general  in  that  it  depends  on  neither  the  spa¬ 
tial  discretization  schemes  nor  the  specific  factorization  form.  Other  implicit 
schemes  than  second  order  backward  Euler  can  also  be  incorporated  readily. 
A  disadvantage  of  this  linearized  scheme  is  the  requirement  of  large  memory 
storage  of  the  L  and  U .  But  when  the  algorithm  is  implemented  on  the  mod¬ 
ern  parallel  computers,  the  demand  for  memory  is  usually  satisfied  without 
major  difficulties. 


3.2  Approximate  Factorization 

As  mentioned  earlier,  approximate  factorization  is  needed  to  reduce  a  multi¬ 
dimension  problem  to  a  sequence  of  one  dimension  problems  in  terms  of 
matrix  inversion.  For  simplicity,  we  consider  the  governing  equations  in  three 
dimensional  Cartesian  coordinates,  where  the  Jacobian  A{U)  in  (3.9)  can  be 
generally  expressed  as 

A{U)  =  AU^  +  BUy  +  CU,  +  DU 

—  VxxU^xx  ~  VxyU^y  —  VxzU^xz 
~  ^yy^.yy  ~  ^vz^,yz  ~  ^zzU^zz  (3-11) 

here  the  matrices  (A,  B,  C,  Z),  Vij)  are  functions  of  U  and  its  gradients,  (see 
for  example  Collis  (1997)).  Substitute  (3.11)  into  the  left  hand  side  of  (3.9) 
and  factorize  it  into  three  spatial  directions,  we  obtain  the  factorized  iteration 


scheme  as  follows 

[5  -h  AriAAx  +  D-  V;^A„)]P  =  - 

(3.12) 

[S  -h  AriBAy  -  VyyAyy)]Q  =  SV 

(3.13) 

[5  -f-  Ar(CA,  -  VxzAxz)]AU’‘  =  SQ 

(3.14) 

where 

2At 

(3.15) 

and  Ax,Axx  ... 

are  the  spatial  differencing  operators. 
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We  can  further  examine  the  eifectiveness  of  the  factorization  by  r  ewriting 
(3.12  -  3.14)  as  a  general  implicit  operator  £(f7") 

£(tr)AC/*  =  -^(31/'=  -  4i7”  +  cr*-^) 


2At' 
-AtF(C/*) 


(3.16) 


where 


C{XP)  =  [5  +  Ar(^A^  +  D  -  14Ax)]5"^ 

X  [S  +  At(SA,  -  n,A^)]5-i 
X  [S  +  At(CA,  -  y„A«)]  (3.17) 

is  a  function  of  f/"  only.  Introducing  the  subiteration  residual  ^  - 

U  where  U  denotes  the  converged  solution  of  subiteration,  e*  satisfies  the 
following  equation: 

£([/") (e*+'  -  e*)  =  -^(3^  - 

3Ar  1. 


2  At 


-At[F{U)  +  A{U)e'‘]  - 


2At 


(3.18) 


Since  U  is  the  converged  solution,  the  above  equation  is  simplified  as 

.A  A  /V7\l 


£{tr)(e‘+'  -  «*)  =  +  i^rA{U)]e 


It  follows  that 


3Ar 


(3.19) 


6*+^  =  C-\ir)[C{Lr*)  -  ^  -  Ar.4(t/)]e* 


(3.20) 


Thus  the  requirement  for  stability  and  convergence  of  the  subiteration  scheme 
is  that 


II  C-\U^)[C{U'*) 

which  implies  that 

£(t/") « 


3Ar 

2At 

3Ar 


-  AtA{U)]  II  <  1 


2At 


+  AtA{U) 


(3.21) 


(3.22) 


Therefore,  in  order  to  achieve  rapid  convergence,  the  factorized  formulation 
C  needs  to  resemble  the  original  unfactorized  form  as  accurately  as  possible, 
i.e.  keeping  the  factorization  error  minimum.  On  the  other  hand,  as  long  as 
the  subiteration  converges,  we  obtain  the  second  order  temporal  accuracy  by 
(3.2);  the  particular  form  of  the  left  hand  matrix  is  only  important  in  the 
sense  that  subiteration  convergence  can  be  achieved. 
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3.3  Spatial  Discretization 

In  all  the  present  computations,  the  spanwise  direction  z  is  assumed  to  be 
homogeneous.  So  to  compute  the  spatial  derivatives,  the  physical  domain 
(x,  y,  z)  is  first  transformed  into  a  computational  domain  (^,  i),  z)  by  a  two 
dimensional  mapping  (x,y)  with  both  ^  and  t)  belong  to  [0  1].  Con¬ 

sider  derivative  operators  in  the  uniform  computational  space  {^,77),  where 
represent  the  nodal  locations  with  =  A^{i  —  1)  for  1  <  i  <  and 
Tjj  —  Ar){j  —  1)  for  1  <  j  <  N^.  we  need  only  to  present  the  differencing 
schemes  in  ^  direction  as  it  is  completely  analogous  in  the  r/  direction. 

At  interior  nodes,  the  fourth  order  central  difference  scheme  is  used  for 
first  and  second  derivatives. 

{%)  =  -  fw)]  (3.23) 

~  ’  2/i  +  /i+2  +  16{/i_i  —  2fi  +  /j+i){p.24) 

Near  the  computational  boundaries,  the  finite  difference  stencil  need  to  be 
biased  toward  the  interior.  As  in  the  interior  points,  a  five  point  stencil  is  used 
at  the  boundary  grid  point  ^i.  The  resulting  difference  schemes  are  fourth  and 
third  order  accurate  for  the  first  and  second  derivatives  respectively.  Hence, 
At  the  first  grid  point  ^1,  the  first  and  second  derivative  are  expressed  as  : 

{%),  "  i^[-25/i  + 48/2  -  36/3  +  16/4  -3/5]  (3.25) 

{w)i  "  +  +  (3.26) 

(3.27) 

while  at  the  second  point  ^2,  they  are 

{%),  "  Y5^|-Vi-10/2+18/3-6/,+/5|  (3.28) 

(0)  =  +  + 

Similar  expressions  hold  for  the  derivatives  at  nodes  of  7V^  —  1  and  but 
with  the  stencils  reversed  and  the  signs  switched  on  the  coefficients  for  the 
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first  derivatives.  With  the  above  fourth  order  discretization  scheme,  the  left 
hand  side  in  each  of  the  (3.12)  -  (3.14)  is  a  block  penta-diagonal  matrix,  see 
the  structures  in  figure  3.2  and  3.3,  with  additional  blocks  at  the  boundary 
nodes  to  account  for  high  order  boundary  treatment.  Each  block  is  a  5  x  5 
matrix  for  the  five  flow  variables.  The  resulting  coeflScient  matrix  can  be 
solved  eflBciently  using  direct  Guassian  elimination  method. 

After  the  derivative  is  obtained  in  the  computational  space,  they  are 
converted  to  the  physical  space  through  the  metrics  of  the  mapping  function. 
For  first  order  derivative,  it  follows 


and  similarly  for  second  derivative 


(3.30) 


^  dx&y 

((  \2 


^,x'0,y  "I"  ^,yV^ 
2  V,y 


'n,x'n,y 

{V,yy 


dri^  ) 


s,xx 

'n,xx 

^,xy 

^,yy 

V,yy  . 

{ 


(3.31) 

where  (^^,  ...)  are  metrics  of  the  mapping  transformation.  Most  of  map¬ 

ping  functions  are  given  in  an  explicit  form  of  x  =  x(^,  rf),  y  —  y(^,  t])  instead 
of  ^  =  ^(x,  y),  T]  =  T){x,y).  In  such  cases,  the  metrics  terms  (^,x,  — )  can 
be  obtained  by  solving  the  linear  system  of  equations  resulting  simply  by 
substituting  the  coordinates  x,  y  into  the  operators  (3.30)  and  (3.31).  Com¬ 
putational  experience  also  shows  that  the  x^,  y,, ...  in  the  metrics  expression 
should  be  computed  using  the  same  differencing  scheme  as  those  for  the  flow 
variables  to  achieve  improved  accuracy  (Fletcher,  1991) 

Since  the  subiteration  is  in  fact  a  steady  state  problem,  to  accelerate  the 
subiteration  convergence,  the  pseudo  time  step  size  At  can  be  allowed  to 
vary  spatially.  In  a  general  curvilinear  coordinates,  taking  into  account  the 
physical  domain  of  dependence  within  a  computation  cell,  the  At  are  chosen 
by 


At  =  CFL 


A^  Ay 
A 


(3.32) 


where 


A  =  iJijttjlAy-l- +  c 

X  [( +  (4  + 


(3.33) 
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At;  are  the  cell  size  in  the  computational  domain;  c  is  the  mean  speed  of 
sound.  This  definition  of  A  includes  the  convective  and  acoustic  phenomena 
but  ignores  viscous  diffusion.  The  CFL  in  (3.32)  is  specified  by  a  uniform 
physical  time  step  At  for  unsteady  computations  which  is  independent  of  At 
and  should  reflect  the  time  resolution  requirements  for  the  physical  process 
of  interest. 

Figure  3.4  shows  a  typical  comparison  of  the  convergence  history  between 
the  present  and  standard  subiteration  schemes.  After  first  subiteration  step, 
the  present  scheme  store  the  LU  decomposition  of  the  block  penta  diagonal 
matrix  and  use  them  in  the  subsequent  steps.  For  a  four  order  of  magnitude 
of  residual  reduction,  a  factor  of  4  —  5  speed-up  is  achieved  in  comparison 
with  the  standard  subiteration  scheme.  Note  the  jump  of  residual  at  last 
points  of  the  two  curves  correspond  to  the  beginning  of  the  next  physical 
step. 


3.4  Boundary  Conditions 

Boundary  conditions  are  introduced  to  replace  the  governing  equations  at 
the  inflow,  outflow  and  the  wall  boundary  of  the  computational  domain,  see 
figure  3.1. 

Consider  an  arbitrary  boundary  constraint  at  time  level  n-l- 1  on  the  flow 
variable  U  =  {p, «,  u,  w,  T) 

=  0  (3.34) 

The  general  implicit  treatment  of  boundary  conditions  in  terms  of  8U  - 
(/"+1  _  (jn  jjg  written  as 

(dB\  ^ 

w)  (3-35) 

For  example,  the  no-slip  and  isothermal  boundary  conditions  are  always  en¬ 
force  at  the  wall,  for  which  (3.35)  becomes 

6u^  =  -<  =  0,  5T^  =  -T^  =  -Tr,  (3.36) 

•  The  density  at  the  wall  is  obtained  by  solving  the  continuity  equation. 

The  boundary  conditions  at  inflow  are  need^  to  provide  both  the  far  field 
mean  flow  and  to  introduce  disturbances  or  free  stream  turbulence  in  to  the 
computational  domain.  According  to  characteristic  analysis,  for  a  subsonic 
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flow,  four  incoming  quantities  must  be  specified  with  one  outgoing  quantities 
computed  from  the  interior  domain.  The  specific  choice  of  these  quantities 
depends  on  the  formulation  of  each  problem.  Here  we  constrain  the  entropy, 
spanwise  and  tangential  velocity,  and  the  incoming  Pliemann  invariant.  The 
outgoing  Riemann  invariant  are  computed  by  first  order  extrapolation  from 
the  neighboring  points  inside  the  boundary.  The  locally  one-dimensional 
Riemann  invariants  are  defined  in  the  direction  normal  to  the  inflow  boundary 
as 

„  2c  „  2c  ,  . 

Rl  =  Vn  — ,  ~  ”1  T  (3.37) 

7-1  7-1 

where  c  is  the  local  speed  of  sound  and  is  the  velocity  normal  to  the 
boundary.  So  at  the  inflow  boundary, 

S  =  Soo,  W  =  W',  Vt=Tt-\- 

Rl  =  Rl  +  Rfi,  R2  =  R2int  (3.38) 

here  the  overline  quantities  denote  the  base  flow  in  the  absence  of  inflow 
disturbance;  primed  quantities  represent  the  disturbance  to  be  introduced 
into  the  domain,  and  the  subscript  t  denotes  the  tangential  velocity  along 
the  inflow  boundary  in  the  x—y  plane.  In  the  present  configuration,  the  inflow 
is  place  close  to  the  body  to  reduce  computational  cost,  and  the  flow  passage 
is  chosen  to  have  significant  blockage  to  mimic  the  wind  tunnel  experiments, 
so  the  potential  solution  should  not  be  used  directly  to  form  the  incoming 
Riemann  invariant  because  of  the  development  of  the  boundary  layer  on  the 
body.  Therefore,  we  first  compute  the  two  dimensional  steady  base  flow  using 
the  same  set  of  the  inflow  conditions  but  in  the  absence  of  disturbance.  The 
velocities  Vt  a^id  are  determined  by 

CJ  =  0,  Vn=  Vnp  (3.39) 

Here  u  is  the  vorticity  at  the  inflow  and  Vnp  denote  the  normal  velocity 
obtained  Jrom  the  potential  solution.  These  overline  quantities  are  used  to 
form  the  Ri  and  the  primed  quantities  are  interpolated  from  the  precomputed 
disturbance  flow  field.  The  interpolation  is  implemented  using  fourth  order 
B-spline.  The  Rzint  are  computed  from  the  interior  domain  by  first  order 
extrapolation  Rant  =  “^R^n-i  —  R2n-2-  When  the  Riemann  invariants  are 
obtained,  the  values  of  Vn  and  c  are  obtained  as  follows. 

^^n  =  2^^^”  ^  ~  "^”2 — (3.40) 
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At  the  outflow,  we  follow  Collis  (1997)  and  use  the  parabolized  Navier- 
Stokes  equations,  i.e.  the  streamwise  second  order  derivatives  in  the  equa¬ 
tions  are  neglected.  In  addition,  the  streamwise  pressure  gradient  is  provided 
by  the  corresponding  potential  solution  at  the  exit.  This  treatment  is  found 
to  yield  a  adequate  and  stable  outflow  boundary  conditions  for  the  laminar 
mean  flow  and  turbulence  computations. 

In  the  spanwise  direction,  periodic  boundary  condition  is  always  imposed 
in  the  present  study. 


3.5  Numerical  damping 

The  function  of  numerical  dissipation,  either  implicitly  contained  in  the  dif¬ 
ferencing  scheme  like  upwind  scheme,  or  explicitly  added  to  the  discretized 
equations,  is  to  damp  the  unresolved  high  wave  number  components  gen¬ 
erated  by  nonlinear  flow  interaction.  Because  of  the  symmetric  stencil,  the 
center  finite  differencing  scheme  in  the  present  spatial  discretization  is  by  it¬ 
self  nondissipative  .  In  large  eddy  simulation  of  turbulence,  these  high  wave 
number  components  are  constantly  generated  and,  if  left  uncontrolled,  leads 
to  to  numerical  instability.  To  overcome  the  problem  and  prevent  the  un¬ 
resolved  components  from  corrupting  the  numerical  solution,  the  following 
numerical  dissipation  terms  are  added  to  the  right  hand  side  of  the  equation: 

D  =  +  At^)  (3.41) 

where  is  the  amplitude  of  the  added  dissipation,  and  0  <  <  1  control 

and  range  of  the  dissipation.  Typically  in  laminar  computation,  such  numer¬ 
ical  dissipation  is  not  needed  because  the  spatial  scales  of  the  flow  is  to  some 
extent  known  a  priori;  grid  can  usually  be  generated  with  enough  resolu¬ 
tion.  But  numerical  dissipation  is  generally  needed  in  turbulence  simulations 
to  ensure  numerical  stability.  But  care  must  be  taken  to  ensure  the  added 
numerical  dissipation  is  minimum  and  it  does  not  deteriorate  the  resolved 
solution.  For  this  purpose,  we  choose  cj  in  such  way  that  the  magnitude  of 
the  added  terms  is  significantly  smaller  than  the  that  of  the  truncation 
error  of  the  underlying  spatial  difference  scheme.  Based  on  modified  wave 
number  analysis,  this  is  achieved  by  choosing 

e  <  O.OlNg  (3.42) 
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where  Ng  =  {N^  +  Nr,  +  Nz)/S  is  the  averaged  number  of  grid  points  in  one 
spatial  direction.  A  more  detailed  discussion  on  this  matter  can  be  found  in 
the  appendix  B. 


3.6  Potential  Flow  solution 

Two  dimensional  potential  solution  is  needed  in  the  present  computation  to 
provide  necessary  initial  and  boundary  conditions  for  corresponding  viscous 
computations.  Let  <f>  be  the  velocity  potential  such  that 

« =  <t>^,  V  =  (3-43) 

the  potential  equation  in  a  nonconservative  form  may  be  written  as 


(1  —  M^)(f>^xx  +  (1  “  ^y)4*,yy  —  0 


(3.44) 


where  =  u/c.  My  =  v/c  and  c  =  y/f/M  with  T  computed  from  the 
isentropic  relation 


T’  1+  -  w2) 


(3.45) 


and  M  =  «oo/cc»  is  the  free  stream  Mach  number.  The  no-penetration 
condition. 

1^  =  0  (3.46) 

on 

At  inflow  and  outflow,  the  potential  function  is  set  to  the  values  correspond¬ 
ing  to  the  free  stream  values.  The  solution  of  the  potential  flow  equation  are 
obtained  using  the  same  approximately  factorized  implicit  method  with  the 
fourth  order  finite  difference  scheme  for  spatial  discretization. 


3.7  Similarity  Solutions 

We  choose  the  self-similarity  analytic  solutions  of  compressible  boundary 
layer  with  heat  transfer  near  a  blunt  leading  edge  (Reshotko  and  Beckwith, 
1957)  to  validate  the  present  numerical  method  in  steady  computations. 
The  detailed  derivations  of  the  self-similar  equations  based  on  Stewartson- 
Illingworth  transformation  are  given  in  the  appendix.  The  boundary  layer 
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velocity  and  enthalpy  profiles  at  a  circular  leading  edge  is  computed  numer¬ 
ically  with  different  Mach  numbers,  Prandtl  number  and  different  temper¬ 
ature  ratio  between  the  wall  the  and  incoming  flow.  used.  As  required  to 
possess  a  similarity  solution,  the  viscosity  is  assumed  to  be  a  linear  function 
of  temperature.  The  computational  grid  in  the  x-y  plane  is  generated  by  us¬ 
ing  an  algebraic  multi-surface  method  (Fletcher,  1991)  which  gives  desirable 
orthogonality  at  the  body  surface  and  at  the  inflow/outflow  boundaries.  The 
grid  is  clustered  towards  the  wall  and  the  leading  edge. 

Figure  3.5  to  figure  3.8  show  the  velocity  and  enthalpy  profiles  of  the 
stagnation  point  boundary  layer  for  Pr  =  1.  The  wall  temperature  is  set 
twice  of  the  total  temperature  of  the  incoming  flow,  i.e.  Tu,/Tq  =  2.  It  can 
be  seen  that  for  both  Mach  numbers,  Ma  =  0.1  and  Ma  =  0.8,  excellent 
agreement  of  the  velocity  and  enthalpy  profiles  are  obtained  between  the 
numerical  and  analytic  solution.  When  Pr  ^  1,  similarity  solution  exists 
in  a  strict  sense  only  when  the  external  velocity  is  zero,  i.e.  Ug  =  0.  This 
corresponds  to  the  stagnation  plane  j/  =  0.  So  figure  3.9  and  3.10  show 
the  enthalpy  profiles  for  both  Mach  numbers  atPr  =  0.7  and  T^/Tq  =  0.5, 
and  again  excellent  agreement  is  reached  between  the  numerical  and  analytic 
solutions. 


3.8  Leading  Edge  Receptivity  to  Sound 

Receptivity  is  defined  as  a  process  by  which  external  flow  disturbances  are 
converted  into  instability  waves  (Morkovin,  1969).  For  flat  plate  boundary 
layer  flow,  the  sound  receptivity  refers  the  generation  of  instability  T-S  waves 
inside  the  boundary  layer  by  ambient  acoustic  waves.  In  low  Mach  number 
case,  this  process  represents  a  change  in  characteristic  flow  length  scale  — 
from  an  acoustic,  fast  travelling,  long-wavelength  mode  in  free  stream  to 
a  convective,  slow-travelling  and  short  wavelength  T-S  wave  mode  in  the 
boundary  layer.  This  conversion  can  be  effectively  realized  when  the  leading 
edge  of  the  flat  plate  is  blunt. 

The  numerical  solution  of  sound  receptivity  in  a  compressible  boundary 
layer  on  a  flat  plate  with  a  super  ellipse  leading  edge  (Lin,  1992;  Collis,  1997) 
is  presented  here.  Unlike  the  usual  numerical  approach  for  receptivity  study 
in  which  the  computation  is  based  on  linearized  governing  equations  about 
a  base  flow,  here  the  full  nonlinear  N-S  equations  are  used  both  for  the  base 
flow  and  the  disturbance  computations.  The  evolution  of  the  disturbance 
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is  obtained  by  subtracting  the  base  solution  from  the  instantaneous  solu¬ 
tion.  By  comparing  with  results  obtained  through  linear  approaches,  this 
computation  provides  a  validation  case  for  the  present  numerical  method  in 
unsteady  flow  simulations. 

The  flow  configuration  and  the  computational  domain,  including  the  sur¬ 
rounding  inflow  and  outflow  sponges,  are  shown  schematically  in  figure  3.11. 
The  inflow  sponge  is  used  to  introduce  acoustic  waves  with  specified  fre¬ 
quency  and  amplitude,  whereas  the  outflow  sponge  damps  spurious  waves 
reflected  from  the  outflow  boundary.  The  geometry  of  the  leading  edge  is 
defined  by  a  super  ellipse: 


+  &  =  1  (3-47) 

a  0 

where  the  a  and  b  are  the  half  major  and  minor  axes  of  the  super-ellipse.  The 
junction  between  the  leading  edge  and  the  flat  plate  is  chosen  to  be  a:  =  0, 
and  the  parameter  m,  controlling  the  smoothness  of  the  junction,  is  taken  as 
m  =  4.  The  reference  length  scale  is  chosen  to  be  a  and  the  aspect  ratio  of 
the  super  ellipse  is  6.  The  stagnation  point  is  therefore  at  a:  =  —  1  and  the 
surface  of  the  plate  is  y  =  6  =  1/6.  Based  on  the  free-stream  velocity  Uoo 
and  o,  the  Reynolds  number  Re  =  Uood/i'  is  14400.  The  free  stream  Mach 
number  is  M  —  0.1.  The  Prandtl  number  is  chosen  as  Pr  =  0.71  and  a 
power  law  /i  =  is  used  for  the  viscosity.  The  wall  is  isothermal  and  the 
temperature  equals  to  the  free  stream  stagnation  temperature.  The  number 
of  grid  points  in  streamwise  direction  is  481  and  in  wall  normal  direction  129. 
The  grid  is  clustered  near  the  wall  and  towards  the  leading  edge  with  about 
35  points  across  the  boundary  layer  and  20  points  per  T-S  wavelength. 

To  obtain  the  correct  disturbance  field,  the  base  flow  must  first  be  a  well 
converged  solution  to  the  N-S  equations.  In  other  words,  the  evolution  of 
the  full  solution  must  be  dominantly  due  to  the  evolution  of  the  disturbance. 
Any  changes  of  the  mean  flow  part  in  the  unsteady  computation  should  be 
negligible  compared  to  the  disturbance  itself.  In  the  present  computation, 
the  mean  flow  is  judged  steady  after  the  residue  has  dropped  8  orders  of 
magnitude  from  its  initial  value.  The  mean  flow  wall  vorticity  a;,,,  and  the 
streamwise  pressure  gradient  dp/ds  at  the  wall  are  shown  in  figure  3.12  and 
3.13,  here  s  is  the  arc  length  of  the  surface  from  the  stagnation  point.  The 
results  are  compared  with  those  from  Collis  (1997)  and  the  agreements  are 
very  good.  In  figure  3.14,  the  streamwise  velocity  profile  at  x  =  2.783  is 
shown  as  a  function  of  the  Blasius  variable  %,  defined  as  7?5  =  yy'«e/(x  +  l). 


33 


Again,  the  present  result  agrees  well  with  Collis  (1997). 

After  the  base  flow  is  obtained,  acoustic  wave  is  introduced  through  the 
inflow  sponge  at  a  frequency  u  =  3.312  and  amplitude  A  =  0.001.  With  the 
mean  flow  Mach  number  M  =  0.1,  the  downstream  acoustic  wave  length  is 
Ao  =  20.833.  In  the  the  sponge  region,  the  following  term  is  added  to  the 
right  hand  side  of  the  N-S  equations, 

S  =  {x, y, z, t)  -  Uref{x, y, z, t) ]  (3.48) 

where  the  C4-e/  is  an  arbitrary  reference  state.  For  the  inflow  sponge  Uj-gf  is 
the  base  flow  plus  the  acoustic  wave,  whereas  for  the  outflow  sponge,  only 
the  base  flow  itself,  a  controls  the  distribution  and  strength  of  the  sponges. 
Following  Collis  (1997),  for  the  inflow  sponge  in  y  direction 

(^{n)  =  A,{- — ^)^>  (3.49) 

Vt- Vs  ' 

where  y,  and  yj  are  the  locations  where  the  sponge  starts  and  ends.  (T(y)  = 
0  for  y  outside  (yi,yi).  The  outflow  sponge  in  ^  direction  is  constructed 
similarly.  In  the  present  computation,  TV,  is  chosen  to  be  3;  A,  is  50  for  the 
inflow  sponge  and  100  for  the  outflow  sponge.  The  time  marching  scheme  is 
the  second  order  backward  Euler  scheme  with  linearized  sub-iterations.  For 
the  unsteady  computation.  At  is  chosen  to  be  around  1000  time  steps  per 
acoustic  period  and  four  sub-iterations  are  used  in  each  time  step. 

The  overall  disturbance  flow  fields  in  the  boundary  layer  contains  not  only 
the  excited  T-S  wave,  but  also  the  incoming  and  scattered  acoustic  waves. 
To  extract  the  T-S  wave,  we  need  first  remove  the  the  acoustic  component 
from  the  total  solution.  In  low  mach  number  flows,  such  a  procedure  was 
introduced  by  Wlezien  (1994).  The  basic  idea  is  that,  in  low  Mach  number 
flows  the  effect  of  the  acoustic  wave  can  be  considered,  due  to  its  long  wave 
length  compared  with  the  T-S  wave,  as  to  change  the  center  of  the  circular 
orbit  in  the  velocity  phase  plane  for  the  T-S  wave.  By  tracking  and  sub¬ 
tracting  this  center  movement,  the  correct  amplitude  of  the  T-S  wave  can  be 
recovered.  This  procedure  works  well  for  the  T-S  waves  on  the  flat  portion 
when  vertical  velocity  is  used  to  denote  the  T-S  wave  signal,  but  does  not 
work  in  the  curved  leading  edge  region,  or  in  high  Mach  number  flows  where 
the  wave  lengthes  for  the  acoustic  and  T-S  waves  are  comparable. 

An  instantaneous  disturbance  field  is  shown  in  figure  3.15  (excluding  the 
sponge  region)  after  the  flow  has  reached  a  periodic  state  driven  by  the  acous¬ 
tic  wave.  The  incident  acoustic  wave  are  clearly  shown  from  p  and  u  contours 
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away  from  the  plate.  The  v  contours  show  the  generation  and  downstream 
development  of  the  T-S  wave  inside  the  boundary  layer.  Generated  at  the 
leading  edge,  the  T-S  wave  pattern  becomes  clearly  visible  on  the  flat  portion 
of  the  wall,  where  the  vertical  component  of  the  acoustic  velocity  is  negligi¬ 
ble.  In  the  leading  edge  region  itself,  however,  the  scattered  acoustic  wave 
is  of  the  same  or  greater  order  of  magnitude  than  the  T-S  wave;  the  two 
velocity  components  can  not  be  cleanly  separated. 

Figure  3.16  shows  the  local  maximum  of  the  amplitude  of  vertical  compo¬ 
nent  of  total  disturbance  velocity  and  the  extracted  acoustic  component  as 
functions  of  x.  The  solutions  from  Collis  (1997)  are  also  plotted  for  compar¬ 
ison.  The  total  disturbance  amplitudes  follow  almost  identical  shape  along 
X,  with  the  present  solution  slightly  lower.  The  acoustic  wave  components 
compares  very  well  throughout  out  the  domain. 

By  subtracting  the  acoustic  components  from  the  total  disturbance  solu¬ 
tion,  the  amplitude  of  the  T-S  wave  is  shown  in  figure  3.17.  The  results  from 
Lin  (1992)  for  incompressible  flow  and  from  Collis  (1997)  are  also  plotted. 
Again,  we  can  see  that  three  computations  follow  similar  shapes  with  some 
small  quantitative  differences.  Given  the  great  sensitivity  of  the  growth  rate 
of  the  T-S  wave,  the  difference  is  considered  to  be  acceptable,  and  the  overall 
agreement  is  satisfactory. 
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Inflow 


Figure  3.1:  Flow  configuration  and  inflow,  outflow  and  wall  boundary 
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Figure  3.2.  The  structure  of  block  penta-diagonal  coefficient  matrix  resulting 
from  implicit  fourth  order  finite  difference  scheme.  Each  block  is  a  5  x  5 
matrix. 
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Figure  3.3:  The  structure  of  periodic  block  penta-diagonal  coefficient  matrix 
resulting  from  implicit  fourth  order  finite  difference  scheme.Each  block  is  a 
5x5  matrix. 
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CPU  time  units 

Figure  3.4:  Comparison  of  the  convergence  history  between  the  present 

and  standard  subiteration  schemes.  - :  present  linearized  subiteration. 

:  standard  subiteration,  symbols  stand  for  each  subiteration  step. 
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Figure  3.6:  Comparison  of  boundary  layer  enthalpy  profile  at  the  leading 

edge.i^e  =  10®,  M  =  0.15,  T^/To=2.0,  Pr  =  1.0. - :  computation, 

- :  similarity  solution 


Figure  3.7:  Comparison  of  boundary  layer  velocity  profile  at  the  leading  edge. 

Re  =  10®,  M  =  0.8,  Tyj/To  =  2.0,  Pr  =  1.0. - :  computation, - : 

similarity  solution 
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Figure  3.8:  Comparison  of  boundary  layer  enthalpy  profile  at  the  leading 

edge.ile  =  10*,  M  =  0.8,  Tu,/To=2.0,  Pr  =  1.0. - :  computation, 

- :  similarity  solution 
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Figure  3.9:  Comparison  of  boundary  layer  enthalpy  profile  at  the  leading 

edge.iJe  =  10*,  M  =  0.15,  T,„/To=0.5,  Pr  =  0.7. - :  computation, 

- :  similarity  solution 
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Figure  3.10:  Comparison  of  boundary  layer  enthalpy  profile  at  the  lead¬ 
ing  edge.ile  =  10^  M  =  0.8,  T,„/ro=0.5,  Pr  =  0.7. - :  computation, 

- :  similarity  solution 
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Figure  3.11:  Schematic  of  computational  domain  and  the  inflow  and  outflow 
sponges. 
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Figure  3.12:  Pressure  gradient  dplds  along  the  plate, 
putation, - ,  result  of  Collis  (1997). 


present  com 


Figure  3.13:  Comparison  of  the  wall  vorticity 
tation, - ,  result  of  Collis  (1997). 


present  compu- 


a) 


b) 


c) 


Figure  3.15:  Disturbance  flow  flelds  of  sound  receptivity  on  a  flat  plate  with  a 
super-ellipse  leading  edge,  a)  density  p,  b)  streamwise  velocity  u,  c)  vertical 
velocity  V.  Contour  levels:  p^ax  =  5.07  x  IQ-^,  =  -9.32  x  Ao  - 

7.59  X  10-.  =  4.80  X  10-,  =  -9.20  x  10-,  An  =  6.55  iTo'^ 

^max  —  117  X  10  Vmin  =  “3.90  X  10“^,  Av  =  8.16  X  10~®. 
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Figure  3.16:  Local  maximum  of  the  amplitude  of  vertical  component  of  dis¬ 
turbance  velocity  and  the  acoustic  components.  - :  total  disturbance. 

- :  acoustic  component. - and . :  results  of  Collis  (1997). 


Figure  3.17:  Local  maximum  of  the  amplitude  of  the  T-S  wave  based  on  the 

vertical  disturbance  velocity.  - :  present  computation. - :  result 

of  Collis  (1997). - :  result  of  Lin  (1992). 
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Chapter  4 

Generation  of  Free  Stream 
Turbulence 


In  this  chapter,  we  describe  the  methodology  of  generating  free  stream  tur¬ 
bulence  for  the  main  large  eddy  simulations.  In  order  to  avoid  the  artificial 
periodicity  of  recycling  the  same  free  stream  turbulence  data,  a  series  of  in¬ 
dependent  but  statistically  identical,  homogeneous  and  isotropic  turbulence 
data  sets  are  first  computed.  With  the  same  specified  length  scale  and  in¬ 
tensity,  they  are  concatenated  by  a  blending  procedure  to  form  a  sufficiently 
long  time  record  for  the  main  LES.  Ihrbulence  statistics  before  and  after  the 
blending  shows  that  the  blending  procedure  is  a  simple  yet  effective  way  of 
generating  free  stream  turbulence. 


4.1  Inflow  Turbulence 

Consider  compressible  turbulent  flow  in  an  arbitrary  domain  ft  (  as  shown  in 
figure  4.1  )  with  a  nominally  uniform  mean  flow  outside  ft.  The  coordinate 
system  (x,  y,  z)  is  such  that  the  mean  stream  is  aligned  with  the  x-axis.  The 
boundary  9ft  of  the  domain  ft  may  be  split  into  three  parts:  a)  an  inflow 
boundary  9ftj,  b)  an  outflow  boundary  9fto,  and  c)  the  wall  boundary 
These  parts  satisfy  the  requirement  9ft  =  9ft,  U9fto(j9ftu,.  Turbulence 
is  convected  into  the  domain  ft  by  the  mean  velocity  U  across  the  inflow 
boundary  9fti. 

The  incoming  stream  is  assumed  to  consist  of  a  “frozen”  turbulent  field 
being  carried  by  the  mean.  The  time  varying  boundary  condition  on  9ftj 
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then  be  obtained  from  the  frozen  field  by  converting,  through  Tay¬ 
lor  s  hypothesis,  the  x-coordinate  to  time  (Lee  et  al.,  1993).  The  problem 
then  reduces  to  accurate  description  of  the  “frozen”  field.  The  simplest  ap¬ 
proach  IS  to  obtain  a  realistic  turbulence  field  by  a  separate  calculation  of 
homogenous  isotropic  turbulence  with  desired  characteristics  (Reynolds  num¬ 
ber,  length  scale,  intensity  etc.).  However,  this  approach  may  turn  out  to 
be  very  expensive.  For  example,  in  the  free-stream  turbulence(FST)-induced 
byp^s  tr^sition  of  a  boundary  layer,  statistical  convergence  requires  at  least 
six  flow  through  times;  the  domain  for  inflow  turbulence  simulation  is  thus 
~  6L,  X  LyX  Lg  in  size.  If  the  resolution  for  the  isotropic  turbulence  is  the 
same  ^  the  mam  simulation,  the  inflow  turbulence  calculation  becomes  six 
times  larger  than  the  problem  of  interest,  making  this  approach  impractical. 

1  his  limitation  may  be  overcome  by  catenating,  in  random  order,  isotropic 
turbulence  in  smaller  boxes  provided  that  the  catenation  is  carried  out  with¬ 
out  significant  modification  of  turbulence  characteristics.  By  exploiting  isotropy 
and  reusing  individual  boxes,  a  long  time  record  can  be  generated  with  a  lim¬ 
ited  number  of  realizations.  Compared  to  the  brute  force  method,  this  pro¬ 
cedure  is  not  only  cheaper,  but  also  capable  of  generating  inflow  turbulence 
with  time  varying  characteristics.  We  first  describe  the  large  eddy  simulation 
of  the  individual  realizations,  and  then  develop  the  the  blending  procedure 
to  concatenate  these  realizations  and  form  the  long  time  record. 


4.2  LES  of  homogeneous  isotropic  turbulence 

The  individual  isotropic  homogeneous  turbulence  data  are  computed  by  sep¬ 
arate  large  eddy  simulations.  The  numerical  LES  code  is  adapted  from  an 
casting  DNS  code  (Lui  and  Lele,  2001)  by  incorporating  the  same  dynamic 
SGS  model  d^cnbed  in  Chapter  2.  The  temporal  and  spatial  dicretization  is 
fourth  order  Rung^Kutta  method  and  sixth  order  compact  finite  difference 
^heme  (Lel^  1992).  Periodic  boundary  conditions  are  applied  in  all  spatial 
directions.  The  initial  condition  is  a  solenoidal  velocity  field  with  uniform 

density  and  temperature  field.  The  initial  three  dimensional  energy  spectrum 
IS  of  the  form 

E{k)  cx  k*  exp[-2(«//Cp)2]  (4  1^ 

where  the  peak  wave  number  Kj,  is  equal  to  4.  After  first  validating  the 
LES  code  against  both  the  DNS  of  Lee  et  al.  (1991)  and  the  experiment  of 
Comte-Bellot  and  Corrsin  (1971),  we  apply  the  simulation  to  a  rectangular 
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box  of  the  size  1.6Z),  1.6D  and  0.4D  in  x,  y  and  z  direction,  where  D  is  the 
diameter  of  the  elliptic  leading  edge  to  be  placed  in  the  domain  of  the  primary 
simulation.  Note  that  the  spanwise  box  size  0.4Z)  is  chosen  to  be  the  same 
as  in  the  primary  simulation.  To  preserve  the  isotropy  of  the  turbulence, 
the  grid  size  is  chosen  to  be  128x128x32,  i.e.  Ax  =  Ay  =  Az.  A  few  small 
wave  numbers  which  represent  the  long  waves  permissible  in  the  x  and  y 
directions  but  not  in  z  are  zeroed  out  initially  to  preserve  the  turbulence 
isotropy.  After  the  turbulence  statistics  such  as  intensity  and  length  scale 
reach  predetermined  values,  snapshots  of  instantaneous  flow  flelds  are  taken 
as  one  turbulence  realization. 

To  compare  directly  with  experimental  measurements,  throughout  this 
report  the  turbulence  integral  scale  L  is  determined  by  a  least-squares  curve¬ 
fitting  to  the  autocorrelation  data  using  an  exponential  function  (Van  Fossen 
et  al.,  1995) 

R(r)  =  (4.2) 

Data  between  0.33  <  R{t)  <  1.0  were  used  for  curve  fitting.  The  exponential 
function  does  not  have  the  correct  limiting  behavior  for  very  small  values  of 
r,  but  the  fit  is  satisfactory  over  the  main  portion  and  the  fitted  R{t)  can  be 
integrated  from  0  to  oo  to  give  the  turbulence  integral  length  scale.  Figure 
4.2  shows  the  initial  and  final  three  dimensional  energy  spectrum,  where  the 
small  scale  turbulence  has  filled  the  high  wave  number  space,  figure  4.3 
shows  the  velocity  skewness  reaches  the  typical  value  —0.4  —  0.5  for  the 
realistic  turbulence  at  Rcl  about  80.  The  time  development  of  turbulence 
kinetic  energy  is  shown  in  figure  4.4.  Notice  that  vP'  is  slightly  larger  than 

and  in  the  time  series,  which  may  be  caused  by  the  different  size 
of  the  computational  domain  in  the  z  direction  compared  to  the  x  and  y 
directions.  Since  this  lack  of  isotropy  is  small,  it  is  not  expected  to  have  any 
major  effect.  As  we  started  with  the  uniform  density  field.  Figure  4.5  shows 
the  time  development  of  the  RMS  value  of  density  fluctuation.  The  time 
development  of  the  dynamic  model  coefficient  C  is  shown  in  figure  4.6. 


4.3  Blending  Procedure 

Before  we  introduce  the  blending  zone,  in  which  the  two  date  set  transition 
smoothly  from  one  data  to  another,  consider  first  two  independent,  statisti¬ 
cally  identical  random  velocity  fields  and  that  are  homogeneous  and 
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isotropic.  Let  u  be  a  linear  combination  of  the  form 

u  =  (4.3) 

where  a  and  P  are  scalar  constants.  If  and  are  velocity  fields  corre¬ 
sponding  to  homogeneous  isotropic  turbulence,  their  mean  is  zero,  as  is  that 
of  the  new  field  u.  Furthermore,  the  two  point  correlation  /L,  (r)  for  the  new 
field  u,  ^ 


W  -  (u,(x)uj(xH-r))  =  ^(x)«j^^(x-l-r)^  + 13^ 

+  afi  ((«S^^(x)uf  (x-l-r)^  + 

reduces  to, 

Hij{T)  =  {q^  +  (4.5) 

by  virtue  of  and  being  independent,  where  the  angled  brackets  (•) 
indicate  an  averaging  (volume  average  suffices  for  the  homogeneous  fields  un¬ 
der  consideration).  To  retain  the  second  order  statistics  of  the  original  fields 
and  an  appropriate  rescaling  of  the  definition  in  (4.3)  is  required. 
This  renormalization  yields 


u  = 


a 


P 


y/a^  +  P^  y/a^  +  p^ 


u 


u  =  cos  9  -f  sin  9 


(2) 

(4.6) 

as 

(4.7) 

The  new  field  so  obtained  retains  the  mean  values  and  second  order  statistics 
of  the  original  fields. 

This  linear  combination  may  be  generalized  by  varying  9,  allowing  smooth 
transition  from  one  field  to  another  over  a  blending  zone  within  which  9  varies 
from  0  to  f  as  illustrated  in  figure  4.7.  With  9  =  9{x)  and  (•>  restricted  to 
averaging  in  y-z  plane,  the  single  point  statistics  are  preserved,  as  are  two 
point  correlations  in  the  y  —  z  plane.  The  two  point  correlation  in  the  x 
direction,  by  (4.7),  can  be  expressed  as 


<  ni(x)uj(x  -I-  r*)  >  =  cos  [  0(x,  r^)  ]  /2.-,(r*)  (4.8) 
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where  0(x,  r)  =  9{x  +  r)  -  0{x),  and  Rij  is  the  original  two  point  correlation 
in  X.  Let  us  consider  the  case  where  r  <  L,  where  L  is  the  correlation  length 
of  the  original  field,  then  in  the  blending  zone,  (j>{x,  r)  can  be  estimated  as 


where  Li,  =  {'K/2)/{d6ldx)max  is  defined  as  a  characteristic  length  of  the 
blending  zone.  PVom  (4.8),  the  difference  of  the  two  point  correlation  in  x 
between  the  blended  and  the  original  field  is  proportional  to  when  ^  <C  1, 
so  we  have 

\  <Ui{x)uj{x  +  rx)  > -Rij{rj:)\  ~ 

which  means  the  longer  the  blending  zone,  the  less  error  will  result  in  the 
two  point  correlation  in  x  direction. 

The  dependence  of  0  on  x  within  the  blending  zone,  however,  also  intro¬ 
duces  an  extra  term  Dg  in  the  dilatation  field: 

V  -  u  =  cos  9  (V  •  +  sin  9  (V  -  +  Ve 

X>e  =  (— sin^u^^^  -I- cos  01  (4-11) 

In  incompressible  flows,  "Dg  violates  mass  conservation,  whereas  in  compress¬ 
ible  flows  it  can  cause  large  artificial  pressure  fluctuations.  This  undesired 
dilatation  may  be  removed  by  a  projection  method  based  on  the  Helmholtz 
decomposition  theorem  for  the  velocity  vector  u 

u  =  VxA  +  Vv?  (4.12) 

where  A  and  ip  are  the  vector  and  scalar  potential  of  u.  The  scalar  poten¬ 
tial  (fe  corresponding  to  the  extra  dilatation  X>e  satisfies  a  Poisson  equation 
obtained  by  taking  the  divergence  of  (4.7), 

VVe  =  2>e  (4.13) 

This  equation  is  solved  with  homogeneous  boundary  conditions  in  x  direction 
and  periodic  in  the  other  two  directions.  Subtracting  V(pe  from  u  yields  the 
the  final  expression  of  the  blended  velocity  u  as 

u  =  u  —  Vipe  (4.14) 
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Two  sample  data  sets  of  homogeneous  isotropic  turbulence  are  used  to 
demonstrate  the  blending  procedure.  The  two  data  sets  may  be  regarded 
as  different  realizations  of  the  same  turbulent  statistical  state,  to  which  the 
above  blending  procedure  is  applied.  In  the  present  case,  the  size  of  the 
blending  zone  is  taken  as  1/8  of  the  box  length.  The  results  for  blending 
over  1/16  of  the  box  length  are  found  to  be  similar. 

The  energy  spectra  of  the  original  and  blended  turbulence  fields  are  shown 
in  figure  4.8.  The  energy  spectrum  of  the  blended  field  is  very  close  to  the 
original.  Because  of  the  increased  length  in  the  x-direction,  a  lower  wave 
number  with  very  small  amount  of  energy  is  present.  The  effect  of  the  pro¬ 
jection  procedure  is  shown  in  figure  4.9,  where  dilatation,  averaged  over  y-z 
plane,  is  plotted  as  a  function  of  x.  Simple  blending  using  equation  (4.7) 
produces  a  dilatation  spike  inside  the  blending  zone  which  is  removed  by 
(4.14).  Notice  that  outside  the  blending  region,  the  flow  field  remains  intact. 
Numerically,  the  velocity  associated  with  the  extra  dilatation  Vv?e  is  at  least 
,  one  order  of  magnitude  smaller  compared  to  u.  The  effect  of  removing  extra 
dilatation  on  the  overall  second  order  statistics  is  therefore  negligible.  The 
velocity  gradient  skewness  averaged  in  y-z  plane  is  plotted  in  figure  4.10  as 
a  function  of  x.  For  clarity,  the  second  half  of  the  first  data  set  and  the  first 
half  of  the  second  data  set  are  shown  together  on  top.  The  blended  data  of 
the  corresponding  range  is  shown  at  the  bottom.  All  three  components  of 
the  skewness  maintain  similar  magnitude  in  the  blending  zone  without  ex¬ 
cessively  large  fluctuations.  The  corresponding  vorticity  distribution,  shown 
in  figure  4.11  before  and  after  blending,  is  also  well  behaved. 

In  this  Chapter,  a  blending  procedure  is  described  for  combining  realiza¬ 
tions  of  homogeneous  isotropic  turbulence  into  a  unified  field  that  serves  as 
a  realistic  representation  of  free-stream  turbulence.  Different  realizations  are 
catenated  by  a  smooth  blending  function,  and  extra  dilatation  is  removed  us¬ 
ing  Helmholtz  vector  decomposition  theorem.  By  construction,  the  combined 
field  preserves  the  turbulence  intensity,  and  the  change  to  other  statistical 
quantities  are  shown  to  be  minimal.  Examples  will  be  given  from  the  LES  of 
free-stream  turbulence  effect  on  leading  edge  heat  transfer  in  the  subsequent 
chapters.  This  simple  yet  effective  method  could  be  useful  in  other  direct 
or  large  eddy  simulations  in  which  effects  of  sustained  free-steam  turbulence 
are  important. 
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Figure  4.1:  Length  of  free-stream  turbulence  field  needed  for  typical  spatial 
numerical  simulations. 


Figure  4.2:  Initial  and  final  three  dimensional  energy  spectrum  for  isotropic 
homogeneous  turbulence  in  a  rectangular  box. 


X 


Blending  zone - 0  <  e{x)  <tx/2 


Figure  4.7:  Blending  function  for  joining  two  turbulence  realizations. 


Figure  4.8:  Energy  spectra  for  the  original  and  blended  data  sets _ • 

data  1. - :  data  2, - :  blended 
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Figure  4.9:  Divergence  of  blended  compressible  turbulence  data  before  and 

after  the  removal  of  extra  dilatation  field.  - :  original, - ,  extra 

dilatation  removed. 


Figure  4.10:  Velocity  gradient  skewness.  The  blending  zone  are  marked 
between  dashed  lines. - , - ,  and . are  S^g  in  x,  y,  z  directions. 
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Chapter  5 

Numerical  Simulation  Results 


In  this  chapter,  we  present  the  numerical  results  of  LES  of  free-stream  tur¬ 
bulence  impinging  upon  an  elliptical  leading  edge  and  the  related  heat  trans¬ 
fer  enhancement.  The  problem  set  up  is  based  on  the  experimental  results 
of  Van  Fossen  et  al.  (1995),  and  the  LES  results  shows  very  good  agree¬ 
ment  with  the  experimental  measurements.  The  interaction  process  between 
the  free  stream  turbulence  and  leading  edge  is  analyzed  through  the  turbu¬ 
lence  statistics.  The  characteristic  vortical  structures  responsible  for  the  heat 
transfer  enhancement  are  identified. 


5.1  Experimental  Setup 

In  the  wind  tunnel  experiments  of  Van  Fossen  et  al.  (1995),  free-stream  turbu¬ 
lence  is  generated  by  different  grids  with  different  length  scale  and  intensity, 
see  figure  5.1.  The  heat  transfer  is  measured  in  the  stagnation  region  of 
four  elliptical  leading  edge  models;  the  ratios  of  major  to  minor  axes  are  1:1, 
1.5:1,  2.25:1,  and  3:1.  The  Reynolds  number  based  on  diameter  of  curva¬ 
ture  at  leading  edge  and  the  mean  streamwise  velocity  ranges  from  37000 
to  228000.  Turbulence  intensity  (measured  in  the  absence  of  the  models)  is 
in  the  range  of  1.1  to  15.9  percent.  The  ratio  of  turbulence  integral  length 
scale  to  the  leading  edge  diameter  range  from  0.05  to  0.30.  In  this  section  we 
describe  the  specification  of  turbulence  intensity  and  length  scale  according 
to  the  experiments.  We  choose  to  use  experimental  data  set  No.  244,  corre¬ 
sponding  the  leading  edge  with  3:1  major  to  minor  axis  ratio,  as  an  example 
to  describe  the  method.  In  our  simulations,  this  corresponds  to  the  case  B 
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in  table  5.2  where  the  flow  conditions  and  computational  parameters  used  in 
the  present  LES  are  summarized.  The  main  simulation  results  that  will  be 
presented  in  this  chapter  are  also  from  case  B. 

First,  the  Reynolds  number  Re  in  the  simulation  is  defined  as 


— 


_  {pU)ooD 

A^oo 


2(1  -  B)  {pUUD 
2-B  u 


=  0.8284  Re„ 


(5.1) 


where  Rcay  is  the  Reynolds  number  used  in  the  experiments,  and  {pU)^^  is 
the  average  between  the  mass  fluxes  at  far  upstream  and  the  location  with 
maximum  model  blockage.  The  blockage  coefficient  B  (Maximum  model 
thickness  /  Wind  TXinnel  Height)  is  0.293  for  the  3:1  model.  The  experimental 
condition  is  corresponding  to  the  data  set  No.  244  for  this  model  and  the 
measured  turbulence  data  are: 


Reav  =  50700  T«  =  0.0352  A^/Z>  =  0.121,  D  =  6.6  cm  (5.2) 

Note  the  values  of  Tu  and  D  are  taken  at  the  corresponding  location  of 
leading  edge  but  in  the  absence  of  the  model.  The  following  experimental 
correlations  are  obtained  for  Tu  (in  percentage)  and  A*  as  a  function  of  the 
streamwise  distance  x  from  the  turbulence  generating  grid: 

T«  =  a(|r  (5-3) 

where  b  is  the  bar  width  of  the  grid.  For  the  grid  used  in  this  group  of 
measurement,  grid  G3,  the  coefficients  in  equation  (5.3)  are 

a  =  149.4,  b  =  0.318  cm,  m  =  -0.830,  I  =  0.264  (5.4) 

Thus,  corresponding  to  the  experimental  value  given  in  equation  (5.2),  we 
can  obtain  that  the  leading  edge  is  located  at 


^LB  =  91.45  b  =  4.4  D  (5.5) 

downstream  the  grid.  The  corresponding  values  for  Tu  and  at  the  inflow 
boundary  which  is  in  current  simulation  located  1.5D  upstream  of  the  leading 
edge  can  then  be  obtained  by  setting 


xiN  =  xle  -  1.5D  =  2.9D  =  60.19  b 


(5.6) 
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Thus  the  Reynolds  number,  turbulence  intensity  Tu  and  streamwise  integral 
scale  A*  are 

Rez)  =  42000,  Tm  =  0.05,  ^=0-1  (5-7) 

Under  the  operation  condition  of  the  wind  tunnel  at  T  «  300K,  the  kinematic 
viscosity  of  air  is  i/  =  1.60  x  10"®m^/s.  For  the  Reynolds  number  Re  = 
42,000,  the  incoming  upstream  velocity  is  U^o  =  10.18m/s,  which  is  roughly 
corresponding  to  Moo  =  0.03.  To  make  the  compressible  code  run  efficiently, 
we  choose  the  upstream  Mach  number  to  be 

Moo  =  0.15  (5.8) 


5.2  Simulation  procedure 


To  represent  the  experimental  conditions  accurately  with  the  current  flow 
configuration,  the  simulations  are  performed  in  a  progressive  way.  We  first 
generate  the  incompressible  potential  flow  around  the  body  with  blockage 
effect  taken  into  account.  For  this  purpose,  consider  a  point  source  placed 
midway  between  two  parallel  planes  with  a  uniform  incoming  stream  Uoo, 
the  blockage  effect  causes  the  flow  field  to  differ,  particularly  in  pressure 
distribution,  from  the  open  flow  case  where  the  walls  are  absent.  Let  a  be 
the  half  height  of  the  wind  tunnel  and  b  the  half  thickness  of  the  plate  (see 
appendix  A  ),  choose  the  stagnation  point  at  the  center  line  ?/  =  0  as  the 
origin  x  =  0,  the  downstream  blockage  can  be  realized  by  placing  the  point 
source  S  of  strength  2nm,  where 


at  location  x  =  c,  where 


ab 

Uoo 

n{a  —  b) 

(5.9) 

=  ^  In  (1  +  ) 

TT  OtU  oo 

(6.10) 

The  resulting  velocity  field  may  be  expressed  as 

^  ^  MJ„  «>th^  2a  ■ 

2(a  -  b)  cos2(ff )  +  sin2(^)  coth"  b 

y  —  ^00  r _ _ 1 

2(a  -  b)  ^sinh"  ^  +  tan2(f )  cosh"  ^ 


(5.11) 
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Although  in  leading  edge  region,  the  separating  streamline  in  (5.11)  repre¬ 
sents  a  modified  Rankin  body  instead  of  the  exact  elliptic  shape  of  the  model, 
it  is  a  good  approximation  to  the  low  Mach  number  flow  far  upstream  and 
downstream.  So  (5.11)  has  been  used  as  the  initial  condition  for  the  com¬ 
pressible  potential  solution.  The  compressible  potential  flow  is  solved  using 
the  procedure  described  in  Chapter  3  and  the  solution  sought  in  two  steps. 
First  we  solve  the  potential  flow  in  a  large  domain  which  covers  half  of  the 
ellipse  downstream  and  extends  vertically  from  the  bottom  to  the  top  wind 
tunnel  wall.  On  such  a  grid,  the  far  field  boundary  conditions  can  be  read¬ 
ily  provided  by  the  above  incompressible  solution.  After  the  compressible 
solution  is  obtained  on  the  large  grid,  a  smaller  grid  which  only  surrounds 
the  leading  edge  region  is  extracted  for  the  subsequent  viscous  computa¬ 
tions.  Accordingly,  the  boundary  conditions  needed  for  viscous  computation 
on  this  smaller,  inner  grid  can  now  be  provided  through  the  whole  solution 
at  the  boundary,  see  figure  5.2.  This  dual  grid  approach  improves  the  res¬ 
olution  near  the  leading  edge  and  reduces  the  computational  cost  for  three 
dimensional  turbulence  simulation. 

As  aforementioned,  four  different  models  with  the  same  leading  edge  di¬ 
ameter  but  different  major  to  minor  axes  ratios,  from  1:1,  1.5:1,  2.25:1  to 

>  were  used  in  the  experiments  to  produce  different  leading  edge  velocity 
gradient.  These  cases  have  been  duplicated  in  the  compressible  potential 
solution  to  test  the  numerical  accuracy.  Figure  5.4  shows  the  comparison  of 
the  velocity  gradients  at  the  leading  from  the  present  compressible  potential 
solution  at  Afoo  =  0.15  with  the  results  obtained  by  a  panel  method(Van 
Fossen  et  al.,  1995).  Note  the  leading  edge  velocity  gradients  are  affected 
significantly  by  the  blockage  effect,  and  the  excellent  agreement  shows  that 
this  important  feature  of  the  experiments  has  been  captured  accurately. 

For  the  viscous  calculation,  the  flow  parameters  are  Rep  =  42, 000,  Moo  = 
0.15,  T^/Tq  =  1.075,  where  Ty,  and  Tq  are  the  wall  and  total  free  stream 
temperature.  The  laminar  solution  for  density  p,  velocity  u  and  v  and  the 
temperature  T  are  shown  in  Figure  5.5  to  5.8.  The  corresponding  profiles 
along  the  stagnation  streamline  (y  =  0)  are  shown  in  figure  5.9  to  5.11.  It 
can  be  seen  that  while  the  p  and  T  are  only  significantly  changed  inside  the 
boundary  layer,  considerable  gradient  caused  by  the  presence  of  the  model 
exists  in  the  mean  streamwise  velocity  even  at  the  inflow  boundary. 

The  non-dimensional  leading-edge  heat  transfer  coefficient,  or  Frossling 
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number,  is  defined  as 


Nu 

Re'l'‘  "  ~KnATRe'l^ 


(5.12) 


where  is  the  nondimensional  wall  normal  temperature  gradient,  and 
and  kq  is  the  thermal  conductivity  evaluated  at  the  wall  temperature  and 
the  free-stream  total  temperature.  Following  Rigby  and  VanFossen  (1992), 
the  characteristic  temperature  difference  AT  in  (5.12)  is  chosen  to  be  the 
difference  between  the  wall  temperature  and  the  adiabatic  wall  temperature, 
i.e.  AT  =  T^-  Tarv.  The  adiabatic  temperature  can  be  approximated  by 


Ta^  =  y/^{To-Tedge)  +  Tedge  (5.13) 


Here  Tedge  is  the  temperature  at  the  edge  of  the  boundary  layer.  Using  the 
isentropic  relation  and  assuming  constant  pressure  across  the  boundary  layer, 
the  Tedge  is  obtained  by 

Trfpe  =  21  (5.14) 

where  Pq  is  the  free  stream  total  pressure.  The  laminar  computation  of  the 
Fr  distribution  along  the  leading  edge  surface  is  compared  with  experimental 
data  obtained  without  the  turbulence  generating  grids  in  figure  5.12  .  Three 
groups  of  data  at  slightly  different  Reynolds  numbers  (based  on  Uav)  are 
included,  and  the  agreement  is  quite  good.  Note  that  in  such  cases,  the 
background  fluctuation  in  the  wind  tunnel  is  Tu  =  0.3%  and  the  length  scale 
is  L/D  =  2.308. 

The  steady  two  dimensional  viscous  solution  is  then  taken  as  the  LES  base 
flow  to  which  free-stream  turbulence,  generated  by  the  method  described  in 
Chapter  4,  is  introduced  through  the  inflow  boundary. 

However,  prior  to  the  LES  corresponding  to  case  B,  a  preliminary  sim¬ 
ulation,  Case  A,  was  preformed  at  a  lower  Reynolds  number  to  test  the 
numerical  method  and  optimize  the  computation  configuration.  Figure  5.13 
shows  the  turbulence  intensity  along  the  the  stagnation  stream  line.  The 
root-mean-square  values  are  obtained  by  averaging  v*  and  w'  in  time  and 
in  the  spanwise  direction.  The  turbulence  is  largely  decaying  until  it  reaches 
a  distance  of  about  D  from  the  leading  edge  where  the  behavior  of  v'  and 
w'  start  to  change.  Notice  that  close  to  the  body  u'  and  w'  are  amplified 
while  v'  continues  to  decay.  In  Van  Fossen’s  experiment,  a  power  law  curve  of 
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Case 

Ma 

Re 

Tu 

L/D 

Domain  (x,  y,  z) 

Grid 

A 

0.15 

10,000 

0.08 

0.10 

7D  X  lOD  X  D 

384  X  192  X  64 

B 

0.15 

42,000 

0.05 

0.10 

3.5D  X  5D  X0.4D 

192  X  144  X  48 

C 

0.15 

43,740 

0.11 

0.06 

3.5D  X  5D  X0.4D 

192  X  144  X  48 

D 

0.60 

42,000 

0.04 

0.06 

3.5D  X  5D  X0.4D 

192  X  144  X  48 

Table  5.1;  Flow  conditions  and  the  parameters  of  numerical  simulations. 


the  form  Tu  ~  x"*  was  used  to  fit  the  decay  of  free-stream  turbulence  in  the 
absence  of  the  model,  here  x  is  the  distance  downstream  of  the  turbulence¬ 
generating  grid.  The  power-law-fitted  curve  is  also  plotted  in  Fig.  5.13  with 
the  same  exponent  m  =  —0.83  used  in  the  experiments.  The  fairly  good 
agreement  indicates  that  the  initial  decay  rate  of  the  free-stream  turbulence 
has  been  captured  correctly  in  the  simulation. 

Figure  5.14  uses  a  logarithmic  scale  on  a  slightly  shifted  x-axis  to  show 
the  same  data  as  in  figure  5.13  with  the  leading  edge  region  amplified  for 
clarity.  The  transformation  of  the  x-axis  used  here  is  x  =  /op(4.02-x).  Note 
X  =  4  corresponds  to  the  location  of  the  leading  edge.  It  reveals  that  the  free- 
stream  turbulence  goes  through  three  different  stages  as  it  impinges  upon  the 
leading  edge.  The  first  is  free  decay  of  the  turbulence,  corresponding  to  the 
distance  from  the  inflow  to  approximately  x  =  3,  as  the  presence  of  the  body 
has  not  been  felt  strongly  by  the  incoming  turbulence.  The  second  stage  is 
an  inviscid  rapid  distortion  process,  approximately  corresponding  to  the  dis¬ 
tance  from  X  =  3  to  X  =  3.95,  where  the  free-stream  turbulence  experiences 
large  straining  by  the  diverging  mean  streamlines.  A  direct,  quantitative 
comparison  between  the  present  results  and  the  compressible  rapid  distor¬ 
tion  theory  (RDT) {Goldstein,  1978)  is  not  easily  obtained  due  to  the  viscous 
effect  and  the  relatively  small  scale  of  turbulence,  but  the  increase  of 
and  and  decrease  of  are  qualitatively  in  agreement  with  the  tem¬ 
poral  prediction  of  RDT  under  plane  strain  (Batchelor  and  Proudman,  1954). 
The  third  stage  occurs  at  a  distance  very  close  to  the  wall  where  the  viscous 
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dissipation  dominates  and  all  the  turbulence  rapidly  reduces  to  zero  at  the 
wall.  Also  notable  is  that  the  spanwise  velocity  continues  to  increase 
significantly  while  and  start  decreasing  due  to  the  presence  of  the 
wall.  This  is  due  to  inviscid  wall-blocking  or  splatting  effect  (Hunt  and  Gra¬ 
ham,  1978;  Perot  and  Moin,  1995),  in  which  the  wall  normal  fiuctuation  is 
suppressed  by  the  no  penetration  condition  of  the  wall  and  has  to  transfer 
its  energy  to  the  other  two  components.  As  wall  blocking  is  an  inviscid  pro¬ 
cess,  it  typically  extends  beyond  the  viscous  dissipation  range  in  the  normal 
direction. 

It  is  based  on  the  results  of  this  simulation  that  a  dual  mesh  methodology 
is  adopted  for  subsequent  simulations  as  shown  in  figure  5.2.  The  allows  the 
LES  to  focus  on  the  rapid  strain  of  the  free-stream  turbulence  and  the  viscous 
interaction  between  the  turbulence  and  wall.  Hence  in  the  LES  of  case  B,  C 
and  D,  the  inflow  boundary  of  the  mesh  is  chosen  to  be  roughly  corresponding 
to  the  end  of  the  free  decaying  process,  in  the  present  case  at  \.bD  upstream 
the  leading  edge.  The  domain  size  in  the  spanwise  direction  is  chosen  to 
be  0.4D.  Figure  5.3  shows  four  streamwise  locations  along  the  leading  edge 
surface  where  profiles  of  various  turbulence  statistics  at  these  location  will 
be  presented. 


5.3  Inflow  turbulence 

Twelve  independent  realizations  of  homogeneous  isotropic  turbulence  with 
intensity  Tu  =  5%  and  integral  length  scale  C/D  =  0.1  are  pre-computed 
using  the  large  eddy  simulation.  The  LES  code  for  infiow  turbulence  genera¬ 
tion  uses  the  same  dynamic  SGS  model  and  is  adapted  from  a  DNS  code  (Lui, 
2003)  which  incorporates  six-order  compact  finite  differencing  and  fourth  or¬ 
der  R-K  time  marching  scheme.  To  fit  the  aspect  ratio  of  the  main  simula¬ 
tion,  the  LES  of  decaying  turbulence  is  performed  in  a  rectangular  box  of 
size  l.GD,  1.6D  and  0.41)  in  x,  y  and  z  direction,  respectively.  The  number 
of  grid  points  is  128  x  128  x  32  to  ensure  the  turbulence  isotropy  and  initial 
velocity  field  is  solenoidal.  After  the  twelve  independent  turbulence  fields  are 
obtained,  they  are  then  lined  up  spatially  and  joined  together  by  applying 
the  blending  procedure  at  the  interface  between  adjacent  fields,  resulting  in 
a  box  of  turbulence  twelve  times  longer  than  each  individual  realization,  but 
having  the  same  characteristics  as  each  one.  This  long  blended  field  then 
serves  as  the  free-stream  turbulence  which  is  convected  into  the  main  com- 
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putational  domain  through  the  inflow  boundary  located  1.5D  upstream  the 
leading  edge. 

Figure  5.15  shows  the  visualization  of  instantaneous  spanwise  velocity 
within  a  slice  of  ^  =  0  and  the  temperature  gradient  at  the  blade  surface 
near  the  stagnation  point  s  =  0.  It  can  be  seen  that  the  isotropic  turbulence 
eddies  are  strongly  stretched  near  the  leading  edge  region.  As  a  result,  the 
vorticity  in  the  free-stream  in  x  direction  is  compressed  and  reduced, 
but  the  Uy  is  greatly  stretched  and  amplified.  After  the  leading  edge,  these 
vortices  wrap  around  the  leading  edge;  the  vortices  originally  in  transverse 
direction  tend  to  form  strong  streamwise  vortices.  The  strong  and  relatively 
thinly  strectched  vortices  impinging  upon  the  leading  edge  and  modifies  the 
thermal  boundary  layer  significantly.  As  a  result,  the  surface  heat  transfer 
distribution  develops  into  thin,  streamwise  elongated  streaky  structures  as 
shown  in  figure  5.15. 

The  instantaneous  temperature  contour,  corresponding  to  the  flow  field  in 
figure  5.15  are  shown  in  figure  5.16  through  a  series  of  wall  normal  sections  at 
different  streamwise  locations.  Clearly,  the  instantaneous  temperature  con¬ 
tours  develops  into  various  mushroom  like  structures  under  the  influence  of 
the  incoming  turbulence  eddies.  A  sequence  of  visualizations  show  that  these 
mushroom  structures  typically  do  not  occupy  a  fixed  spanwise  location,  but 
moves  laterally  over  a  significant  distance  during  their  lifetime.  The  charac¬ 
terization  of  this  movement  and  its  implication  are  discussed  subsequently. 


5.4  Turbulent  mean  flow  and  heat  transfer 

The  turbulent  mean  density  p,  streamwise  velocity  U,  transverse  velocity 
V  and  temperature  T  are  shown  in  figure  5.17,  5.18,  5.19  and  5.20.  The 
contours  of  the  turbulent  mean  fields  show  qualitatively  similar  patterns 
to  the  corresponding  laminar  quantities,  particularly  outside  the  boundarv 
layer. 

The  profiles  of  turbulent  mean  streamwise  velocity  U  (in  local  s  —  n 
coordinates)  and  temperature  (r-l)/(T,„-l)  are  compared  with  the  laminar 
profiles  at  s  =  0  in  figure  5.21,  0.2I>  in  figure  5.22,  0.8D  in  figure  5.23  and 
1.6Z)  in  figure  5.24.  Of  all  the  locations,  the  turbulent  mean  temperature 
profile  has  a  greater  slope  at  the  wall  while  the  mean  velocity  profile  has 
almost  identical  slope  at  the  wall.  This  means  that  the  heat  transfer  is  much 
more  sensitive  to  the  free-stream  turbulence  than  the  skin  friction. 
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The  time  history  of  spanwise  averaged  FVossling  number  at  the  stagna¬ 
tion  point  is  shown  in  figure  5.25.  Notice  that  the  significant  increase  in  Fr 
around  t  =  4  is  the  time  when  the  free-stream  turbulence  reaches  the  leading 
edge  after  entering  the  domain  at  t  =  0.  The  turbulent  mean  Fr  after  t  =  4 
shows  a  clear  increase  over  the  laminar  value,  indicating  the  heat  transfer  en¬ 
hancement  by  the  free-stream  turbulence.  The  distribution  of  the  Fr  on  the 
surface  is  plotted  in  figure  5.26  along  with  the  experimental  measurements  of 
Van  Fossen  et  al.  (1995)  for  comparison.  Both  the  computation  and  exper¬ 
imental  results  show  that  while  the  surface  distribution  of  the  heat  transfer 
coefficient  is  essentially  the  same  as  in  the  laminar  case,  its  amplitude  is 
significantly  increased  by  the  free-stream  turbulence.  Although  the  present 
LES  result  is  slightly  lower  than  the  experimental  measurements,  as  not  all 
the  turbulence  scales  are  resolved  in  the  computation,  the  overall  agreement 
is  quite  good.  It  will  be  shown  in  Chapter  6  that  except  at  very  small  scales, 
as  the  length  scale  of  the  turbulence  decreases,  its  effect  on  heat  transfer 
increases.  So  it  is  generally  important  to  resolve  the  small  scale  turbulence 
for  good  numerical  prediction  of  stagnation  point  heat  transfer. 

To  quantify  the  effect  of  SGS  modelling,  a  typical  distribution  of  spanwise 
averaged  SGS  i^r  is  shown  in  figure  5.27.  The  distribution  of  time  averaged 
SGS  ut  is  shown  in  5.28.  Both  figures  show  that  the  maximum  ut  occurs 
in  the  leading  edge  region,  where  strong  vortex  stretching  and  amplification 
produces  small  scale  turbulence.  However,  the  value  of  iat  is  relatively  small, 
about  half  of  the  molecular  viscosity,  indicating  the  grid  resolution  close  to 
the  wall  is  adequate. 


5.5  Reynolds  stress  and  turbulent  transport 

Of  main  interest  in  stagnation  point  turbulent  flow  is  the  distribution  of  tur¬ 
bulence  statistics  along  the  stagnation  streamline.  Figure  5.29  shows  stag¬ 
nation  line  distribution  of  RMS  value  of  turbulence  intensity.  As  mentioned 
before,  the  choice  of  the  closer  inflow  boundary  essentially  eliminates  the  free 
decaying  stage  of  the  turbulence.  After  a  relatively  balanced,  isotropic  de¬ 
velopment  stage  before  x  ~  2.4,  the  free-stream  turbulence  becomes  strongly 
anisotropic  due  to  the  mean  flow  straining  effect,  i.e.  Urmsi  Wrms  increase  but 
Vrms  decreases  between  x  2.4  to  x  ~  2.7.  At  this  stage,  the  turbulence 
experiences  the  change  of  mean  flow  due  to  the  presence  of  the  leading  edge, 
but  it  was  not  until  within  even  a  closer  distance  x  >  2.7  that  the  leading 
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edge  surface  itself  exerts  an  strong  blocking  effect  upon  the  turbulence.  The 
blocking  effect  causes  a  continuing,  rapid  increase  in  Wrms  while  Urms  already 
starts  to  decay.  At  the  same  time,  because  the  high  value  of  Umu  and  Wrmg  in 
this  case,  the  previously  decreasing  Vrm$  also  start  to  increasing  after  x  ~  2.7 
in  an  almost  linear  fashion.  Within  an  extremely  closer  distance,  viscous 
effect  dominates  and  all  the  turbulence  vanishes  at  the  wall. 

To  understand  the  dynamics  of  turbulence  Reynolds  stress  and  heat  trans¬ 
port,  it  is  instructive  to  write  their  governing  equations  —  the  Reynolds  av¬ 
eraged  Navier-Stokes  equations.  Let  /  denote  the  time  average  value  of  /, 
we  can  decompose  turbulence  fields  into  mean  and  the  fluctuation 

p  =  p  +  f/,  Ui^Ui  +  u\,  p=P  +  p',  T  =  T  +  r  (5.15) 
and  similarly  for  the  transport  coefficients 

At- =  A -I- A',  Pt  =  P  +  p\  kt  =  k  +  k,',  (5.16) 

Then  the  steady  transport  equation  for  the  Reynolds  stress  is 
pUkiu^u’j)^  = 

-  Ui,k  [p{UkU’j)  +  Ukf/u'j  +  pfu'^u'j] 

-  Uj,k  +  Ukffu'i  +  pfu'f.u\  ] 

-  P<K«i),fc  -  pfu'f^{u\u])^k 

+  ^R  (5.17) 

On  the  right  hand  side  of  (5.17),  the  first  row  is  the  velocity  pressure  gradient 
correlation;  the  second  and  the  third  rows  are  the  Reynolds  stress  production; 
the  fourth  row  is  the  turbulence  transport,  in  the  last  row  is  the  turbulent 
dissipation 

+ A,A^-hA^^)-h 
(ADj^-f  D^  AX  +  A^)] 

t  +  ^ik,k  + 

{P^jkyk^i  ■f’  (5.18) 
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Where  D  =  Ui^i  denotes  the  divergence  of  the  velocity  field.  Similarly,  for 
the  turbulence  heat  flux  u^T',  we  have 

mW%  =  -W^ 

-  U,,t[pu'^T  +  Ut7P  +  P'KT’] 

-  KKrOjt  -  Utp'KT'h  - 

+  $,  +  «!,  (5.19) 

On  the  right  hand  side  of  (5.19),  the  first  row  is  the  temperature  pressure  gra¬ 
dient  correlation;  the  second  and  the  third  rows  are  the  heat  flux  production; 
the  fourth  row  is  the  turbulence  transport  and  the  fifth  row  is  the  turbulent 
heat  conduction.  In  the  last  row  $p,  and  are  the  pressure  dilatation  work 
and  viscous  dissipation.  The  expressions  for  $p  and  are  omitted  here; 
they  are  usually  not  importance  in  understanding  the  turbulence  dynamics. 

Figure  5.30  to  5.32  show  the  turbulence  budget  for  u'^,  v'^  and  w'^  along 
the  stagnation  line.  For  u'^,  there  is  a  large  production  term  corresponding 
to  —u'^dU/dx  due  to  the  mean  flow  straining.  Conversely,  the  production  for 
is  vl^dU/dx  which  has  an  opposite  sign,  and  thus  actually  takes  energy 
away  from  the  vertical  fluctuation.  Since  the  turbulent  mean  flow  is  two 
dimensional,  there  is  not  production  term  for 

The  turbulent  transport  has  the  largest  absolute  value  but  is  mostly  con¬ 
fined  only  within  near  wall  region.  It  also  reverse  sign  as  the  wall  is  ap¬ 
proached.  This  is  characteristic  of  turbulent  transport,  i.e.  when  integrated 
along  the  stagnation  line,  its  net  contribution  to  turbulence  energy  is  zero. 
Turbulent  transport  is  not  important  for  compared  to  other  terms.  For 
turbulent  transport  largely  cancels  with  the  mean  convection  term,  and 
the  summation  of  these  two  plotted  in  5.32  shows  it  has  a  weak  dissipative 
effect  on  the  w'^. 

The  viscous  dissipation  term  extends  from  the  wall  the  largest  distance 
for  u'^  and  smallest  distance  for  It  has  the  largest  value  w'^  yet  smallest 
value  for  u'^.  Except  for  the  dissipation  term  do  not  play  a  significant  role 
at  the  along  the  stagnation  line  unless  within  the  extremely  close  distance  to 
the  wall  for  v'^. 

Of  particular  interest  is  the  redistribution  term,  i.e.  the  correlation  be¬ 
tween  velocity  and  pressure  gradient.  It  has  comparable  amplitude  in  all 
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three  Reynolds  normal  stress  and  plays  very  a  important  role  for  each  of 
them  .  Essentially,  its  overall  effect  is  to  induce  large  turbulence  anisotropy 
in  strongly  strained  mean  flow.  In  figure  5.30,  it  appears  largely  negative  for 
vT-  but  positive  for  xP-  and  shownj^n  figure_5.32  and  5.32.  So  it  takes 

energy  from  vP  and  redistribute  it  to  vP  and  For  both  of  them,  this 
appears  to  be  the  only  major  energy  input.  For  v’^  it  balances  with  the  neg¬ 
ative  production  term,  while  for  it  balances  with  the  viscous  dissipation 
and  total  convection.  It  should  be  pointed  out  that  this  particular  feature  of 
redistribution  term  is  important  for  Reynolds  stress  modelling.  The  turbu¬ 
lence  anisotropy  induced  by  the  redistribution  is  often  the  cause  of  the  failure 
of  using  two  equation  models  in  strongly  strained  flows,  since  in  these  models 
the  turbulence  is  usually  assumed  isotropic  and  only  the  kinetic  energy  k  is 
modelled. 

To  more  clearly  show  the  vortex  stretching  effect,  the  fluctuation  vorticity 
along  the  stagnation  line  is  plotted  in  figure  5.33.  It  can  be  seen  that  the 
amplitude  of  turbulence  vorticity  w'  has  increased  by  two  order  of  magnitude 
from  the  inflow  to  the  wall. 

To  better  visualize  the  typical  vortex  structures  and  its  effect  on  heat 
transfer,  the  temperature  contour  and  the  corresponding  velocity  field  in  the 
stagnation  plane  y  =  0  is  shown  in  figure  5.34  at  a  time  interval  of  At  =  0.6. 
The  velocity  fields  show  clearly  the  formation  of  these  mushroom  shapes  in 
the  thermal  contour;  the  strong,  amplified  y— oriented  vortices  creates  re¬ 
verse  flow  (u  <  0)  in  the  stagnation  region,  which  lifts  up  the  hot  fluid 
particle  close  to  the  wall  and  swaps  them  with  those  away  from  the  wall  at 
a  lower  temperature.  In  doing  so,  the  mushroom  structure  for  temperature 
contour  forms.  However,  it  should  be  noted  that  directly  underneath  these 
mushrooms  structures,  the  thermal  boundary  layer  is  thicker  than  the  undis¬ 
turbed  case  the  heat  transfer  is  in  fact  reduced.  It  is  the  region  between 
these  mushrooms,  where  the  disturbed  flow  has  a  positive  normal  velocity 
towards  the  wall,  the  boundary  layer  become  thinner,  and  consequently  the 
heat  transfer  is  increased.  The  overall  or  spanwise  averaged  Fr  number  ap¬ 
parently  will  depend  on  the  distribution  and  intensity  of  these  thickened  and 
thinned  regions. 

From  figure  5.34,  it  is  observed  that  these  mushrooms  are  moving  along 
the  spanwise  direction.  To  quantify  its  average  convection  speed,  the  space- 
time  correlation  of  the  wall  vorticity  along  the  stagnation  line  (z-axis)  is 
shown  in  figure  5.35.  By  following  the  slowest  descent  line  on  the  contour, 
the  expected  value  of  the  moving  speed  is  found  at  ±0.06,  the  same  order  as 
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Wrms-  The  lateral  movement  of  the  wall  vorticity  can  also  be  explained  by  the 
wall  blocking  effect.  Considering  a  single  vortex  approaching  a  wall  parallel 
to  its  axis,  an  inviscid  image  vortex  with  opposite  sign  will  be  generated 
at  a  equal  distance  on  the  other  side  of  the  wall  in  order  to  satisfy  the  no 
penetration  boundary  condition  on  the  wall.  Then  the  mutual  induction  of 
these  two  vortices  generates  a  velocity  along  the  wall  whose  magnitude  is 
in  general  proportional  to  the  intensity  of  the  approaching  vortex,  and  this 
cause  the  vortex  moves  laterally. 

The  probability  density  function  of  Fr  at  the  leading  edge  is  shown  in 
figure  5.36.  It  deviates  from  a  normal  distribution  and  indicates  that  events 
associated  with  large  positive  heat  transfer  are  less  frequent  than  those  with 
large  negative  variance.  _ _  _ 

Figure  5.37  to  figure  5.46  show  the  distribution  and  profiles  of  //*, 
w'^  and  T'^.  Figure  5.47  to  figure  5.52  show  the  distribution  and  profiles 
of  u'v',  u'w'i  and  v'w'.  Figure  5.53  to  figure  5.58  show  the  distribution  and 
profiles  of  u'T',  v'T'  and  v/T'. 

The  rest  of  the  figures  show  the  magnitude  and  distribution  of  various 
terms  in  the  Reynolds  stress  transport  equation,  including  the  production, 
pressure  strain  and  turbulent  transport.  These  terms  are  typically  important 
for  turbulence  modelling  and  the  profiles  at  specific  locations  are  intended 
to  provide  quantitative  information.  For  brevity,  the  description  of  the  char¬ 
acteristics  of  each  figures  is  omitted  here,  but  the  specific  information  are 
mostly  provided  in  the  captions. 


5.6  High  intensity  and  High  Mach  Number 
cases 

At  the  time  of  this  writing,  the  post  processing  of  the  data  from  the  LES 
of  high  turbulence  intensity,  case  C,  and  high  Ma  number,  case  D,  is  still 
in  progress.  Only  preliminary  results  will  be  presented  here.  Figure  5.83 
shows  the  stagnation  streamline  distribution  of  turbulence  intensity.  It  has 
essentially  similar  features  as  lower  Tu  case  B  except  the  turbulence  decays 
after  entering  the  domain.  This  is  mainly  because  the  smaller  length  scale 
for  case  C.  The  heat  transfer  coefficient  on  the  surface  is  shown  in  figure 
5.84.  It  shows  good  agreement  with  the  experimental  results  but  the  simu¬ 
lation  results  is  also  lower  than  the  measured  data.  Figure  5.85  shows  the 
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Fr  distribution  on  the  surface  for  high  Mach  number  case.  Since  in  this 
case  the  intensity  is  relatively  low,  the  heat  transfer  enhancement  is  modest. 
Qualitatively,  the  shape  appears  similar  to  the  low  Mach  number  cases. 


Figure  5.1:  Wind  tunnel  measurements  for  leading  edge  heat  transfer 
hancement  under  free  stream  turbulence.  (Van  Fossen  et  al.,  1995) 


en- 
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Figure  5.5:  Contours  of  density  p  of  the  laminar  flow  at  Ma  =  0.15,  Re 
42, 000,  T,,/To  =  1.075. 
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Figure  5.7:  Contours  of  transverse  velocity  v  of  the  laminar  flow  at  Ma  = 
0.15,  Re  =  42,000,  T^/Tq  =  1.075. 
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Figure  5.12:  Frossling  number  distribution  on  the  surface.  - :  laminar 

computation,  Symbols:  experimental  measurements  aXTu  —  0.3%,  L/D  = 
2.308;  o:  Re^v  =  50100,  a:  Rbov  —  50200,  □:  Reav  =  50700. 


Figure  5.13:  Turbulence  intensity  along  the  stagnation  streamline. - : 

Wrm», . :  Vrms, - :  lUrms,  - :  power  law  fit  by  a;-®®®.  (Van 

Fossen  et  al.,  1995). 
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Figure  5.14:  Turbulence  intensity  along  the  stagnation  streamline  in  trans¬ 
formed  coordinates  x  for  clarity,  where  x  =  log  (4.02  —  x).  - :  Urms, 
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Figure  5.17:  Mean  density  contours.  Re  =  42, 000,  Ma  =  0.15,  Tu  =  0.05, 
L/D  =  Q.\.  Contour  minimum:  0.9076,  maximum:  1.004,  increment:  0.0048. 
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Figure  5.18:  Mean  streamwise  velocity  contours.  Re  =  42,000,  Ma  =  0.15 

Tu  =  0.05,  L/D  =  0.1.  Contour  minimum:  0.0  ,  maximum:  1.233,  incre^ 
ment:  0.0615. 
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Figure  5.19:  Mean  vertical  velocity  contours.  i?e  =  42,000,  Ma  =  0.15, 
Tu  =  0.05,  LjD  =  0.1.  Contour  minimum:  -0.5926,  maximum:  0.5977, 
increment:  0.0567. 
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Figure  5.20:  Mean  temperature  contours.  Re  =  42,000,  Ma  =  0.15,  Tu  = 
0.05,  L/D  =  0.1.  Contour  minimum  0.995,  maximum  1.075,  increment 


Figure  5.21:  Comparison  between  turbulent  mean  and  laminar  temperature 

profile  along  stagnation  streamline  s  =  0.  - :  turbulent  mean, - : 

laminar. 


Figure  5.22:  Comparison  between  profiles  of  turbulent  mean  velocity  U,  tem¬ 
perature  (T—l)/(Tu,—l)  at  s  =  0.2D  with  laminar  profiles, - :  turbulent 

mean, - :  laminar. 
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Figure  5.23: 
perature  (T- 
mean, - 


Figure  5.24: 
perature  (T- 
mean, - 


Comparison  between  profiles  of  turbulent  mean  velocity  U,  tem- 

-l)/(T,j,— 1)  at  s  =  0.8Z)  with  laminar  profiles, - •:  turbulent 

-:  laminar. 


Comparison  between  profiles  of  turbulent  mean  velocity  U,  tem- 

l)/(7’w-l)  at  s  =  1.6D  with  laminar  profiles, - :  turbulent 

■:  laminar. 
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Figure  5.25:  Time  history  of  spanwise  averaged  Frossling  number  at  the 
leading  edge. 


sfD 

Figure  5.26:  Averaged  Frossling  number  distribution  on  the  body  surface. 

- :  turbulent  mean, - :  laminar,  o  :  experimental  data  (Van  Fossen 

et  al.,  1995). 
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Figure  5.27:  Instantaneous  SGS  eddy  viscosity  vxlv  (averaged  in  z)  .  Re 
42,000,  Ma  =  0.15,  Tu  =  0.05,  LjD  =  0.1.  Contour  minimum  :  O.C 
maximum:  0.75,  increment:  0.014. 
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Figure  5.29:  IXirbulence  intensity  along  the  stagnation  streamline. 

'^rmsy . •  ^rm5>  ~  •  '^rms 


Figure  5.30:  Turbulence  budget  of  along  the  stagnation  streamline. 

°  •  production,  :  dissipation, . :  mean  convection, - : 

turbulent  transport, - :  pressure  gradient  velocity  correlation. 
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Figure  5.31:  Turbulence  budget  of  along  the  stagnation  streamline. 

— o — :  production, - :  dissipation, . :  mean  convection, - : 

turbulent  transport, - :  pressure  gradient  velocity  correlation. 
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Figure  5.32:  Turbulence  budget  of  w'^  along  the  stagnation  streamline. 

- :  dissipation, . :  mean  convection  +  turbulent  transport, - : 

pressure  gradient  velocity  correlation 
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Figure  5.35:  Space-time  correlation  of  vorticity  at  the  leading  edge  (s  =  0). 
Contour  minimum:  0.1,  maximum:  1.0,  increment:  0.002.  The  dash  lines 
indicates  that  the  most  probable  convection  speed  in  spanwise  direction  is 
~  iO.06. 


Fr 

Figure  5.36:  Probability  density  function  of  Frossling  number  at  the  leading 


Figure  5.39:  Turbulent  velocity  fluctuation  Re  =  42, 000,  Ma  =  0.15, 
Tu  —  0.05,  LfD  =  0.1.  Contour  minimum  :  0.0035,  maximum:  0.0233, 
increment:  0.001. 


Figure  5.40:  The  profiles  of  turbulent  velocity  fluctuation  along  wall 

normal  direction. - :  s  =  0, - :  s  =  0.2D, - :  s  =  0.8D,  and 

. :  s  =  1.6D. 
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Figure  5.42:  The  profiles  of  turbulent  velocity  fluctuation  along  wall  nor¬ 
mal  direction.  - :  s  =  0, - :  s  =  0.2D, - :  s  =  0.8D,  and 

. :  s  =  1.6i). 
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Figure  5.43:  TYirbulent  velocity  fluctuation  w'^,  Re  =  42, 000,  Ma  =  0.15, 
Tu  =  0.05,  L/D  =  0.1.  Contour  minimum  :  0.0040,  maximum:  0.0140, 
increment:  0.001. 


Figure  5.44:  The  profiles  of  turbulent  velocity  fluctuation  w‘ 

normal  direction. - :  s  =  0, - :  s  =  0.2D _ s 

. :s  =  1.6D. 


along  wall 
0.8D,  and 
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Figure  5.45:  Turbulent  temperature  fluctuation  Re  =  42,000,  Ma  = 
0.15,  Tu  =  0.05,  L/D  =  0.1.  Contour  minimum  :  0.00002,  maximum 
0.0002,  increment:  0.00002. 
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Figure  5.46:  The  profiles  of  turbulent  temperature  fluctuation  along  wall 

normal  direction. - :  5  =  0, . :  s  =  0.2D, - :  s  =  0.8D,  and 

. :  s  =  1.6D. 


106 


u'v'  X  10^ 


Figure  5.48:  The  profiles  of  turbulent  Reynolds  stress  u'v'  along  wall  normal 

direction.  - :  s  =  0, - :  5  =  0.2D, - :  5  =  0.8D,  and . : 

s  =  1.6D. 


Figure  5.49:  Turbulent  Reynolds  stress  u'w',  Re  =  42,000,  Ma  =  0.15. 
Tu  =  0.05,  L/D  =  0.1.  Contour  minimum  :  -0.0005,  maximum:  0.0005. 
increment:  0.00005. 


Figure  5.50: 
direction.  — 
s  =  1.6D. 
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The  profiles  of  turbulent  Reynolds  stress  v!w'  along  wall  normal 
- :  5  =  0, - :  5  =  0.2D, - :  s  =  0.8Z),  and . : 
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Figure  5.51:  Turbulent  Reynolds  stress  v'v/,  Re  =  42,000,  Ma  —  0.15. 
Tu  =  0.05,  L/D  —  0.1.  Contour  minimum  :  -0.0006,  maximum:  0.0006 
increment:  0.00005. 
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Figure  5.52:  The  profiles  of  turbulent  Reynolds  stress  v'w'  along  wall  normal 
direction.  - :  s  =  0, - ;  s  =  0.2D, - :  s  =  0.8D,  and . : 
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Figure  5.54:  The  profiles  of  turbulent  heat  flux  u’T'  along  wall  normal  di¬ 
rection.  - :  s  =  0, - :  s  =  0.2D, - :  5  =  0.8D,  and . : 

s  =  1.6D. 
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Figure  5.58:  Profiles  of  turbulent  heat  flux  w'T'  along  wall  normal  direction. 
- :  5  =  0, - :  5  =  0.2D, - :  5  =  0.8D,  and . :  s  =  1.6D. 
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Figure  5.60:  Profiles  of  turbulence  production  for  u'^  along  wall  normal  di¬ 
rection.  - :  5  =  0, - :  5  =  0.2D, - :  s  =  0.8D,  and . : 

s  =  1.6D. 
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Figure  5.62:  Profiles  of  turbulence  production  for  elong  wall  normal  di¬ 
rection.  - :  5  =  0, - :  s  =  0.2D, - :  s  =  0.8D,  and . : 

s  =  1.6D. 
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Figure  5.63:  Turbulence  production  for  u'T',  Re  =  42,000,  Ma  =  0.15, 
Tu  =  0.05,  L/D  =  0.1.  Contour  minimum  :  -0.007,  maximum:  -0.0002, 
increment:  0.0002. 


Figure  5.64:  Profiles  of  turbulence  production  for  u'T  along  wall  normal 

direction.  - :  s  =  0, - :  s  =  0.2D, - :  s  =  0.8D,  and . : 

s  =  1.6D. 
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Figure  5.66:  Profiles  of  turbulence  production  for  v'T  along  wall  normal 

direction.  - :  s  =  0, - :  s  =  0.2D, - :  5  =  0.8D,  and . : 

s  =  1.6D. 
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Figure  5.67:  Pressure-strain  correlation  Re  =  42,000,  Ma  —  0.15 

Tu  =  0.05,  L/D  =  0.1.  Contour  minimum  :  -0.115,  maximum:  0.08,  incre¬ 
ment:  0.005. 


Figure  5.68:  Profiles  of 

direction.  - s  =  0 

s  =  1.6D. 


Figure  5.69:  Pressure-strain  correlation  =  42,000,  Ma  =  0.15 

Tu  =  0.05,  L/D  —  0.1.  Contour  minimum  :  -0.08,  maximum:  0.115,  incre¬ 
ment:  0.005. 
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Figure  5.70;  Profiles  of  pressure  strain  correlation  j/S22  along  wall  normal 

direction. - :  s  =  0, - :  $  =  0.2D, - ;  s  =  0.8D,  and . : 

s  =  1.6Z>. 


Figure  5.71:  Pressure-strain  correlation  Re  =  42,000,  Ma  =  0.15, 

Tu  =  0.05,  LfD  =  0.1.  Contour  minimum  :  -0.016,  maximum:  0.023, 
increment:  0.001. 


Figure  5.72:  Profiles  of  pressure  strain  correlation  j/S'^^  along  wall  normal 

direction.  - s  =  0, - :  s  =  0.2D, - :  s  =  0.8D,  and . : 

s  =  1.6£>. 


Figure  5.73:  Turbulent  transport  for  u'^.  Re  =  42,000,  Ma  =  0.15,  Tu  = 
0.05,  L/D  —  0.1.  Contour  minimum  :  -0.05,  maximum:  0.04,  increment: 
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Figure  5.77:  Turbulent  transport  for  v/^,  Re  =  42,000,  Ma  =  0.15,  Tu  — 
0.05,  L/D  =  0.1.  Contour  minimum  :  -0.028,  maximum:  0.012,  increment: 
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2d{jpfu’)/dx,  Re=  42,000,  Ma  =  0.15, 
tiimum  :  -0.16,  maximum:  0.24,  incre- 


Figure  5.80:  Profiles  of  pressure  velocity  correlation  2d{j[/u')ldx  along  wall 

normal  direction. - :  s  =  0, - :  s  =  0.2D, - :  s  =  0.8D,  and 

. :  s  =  1.6D. 
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Figure  5.82:  Profiles  of  pressure  velocity  correlation  2d{p'v')/dy  along  wall 

normal  direction. - :  s  =  0, - :  s  =  0.2D, - :  s  =  0.8D,  and 

. :  s  =  1.6D. 


Figure  5.83:  Turbulence  intensity  along  the  stagnation  streamline  Tu  = 
0.11,  L/D  =  0.06,  Re  =  43740, . :  Ur 
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Figure  5.84:  Averaged  Prossling  number  distribution  on  the  body  surface. 

- :  turbulent  mean, - :  laminar,  o  :  experimental  data  (Van  Fossen 

et  al.,  1995). 
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Figure  5.85:  Averaged  Frossling  number  distribution  on  the  body  surface  at 
Ma  =  0.6. - :  turbulent  mean, - :  laminar 
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Chapter  6 

Theoretical  Analysis 


6.1  Introduction 

The  study  of  fluid  flow  and  heat  transfer  at  a  perturbed  two-dimensional 
forward  stagnation  point  provides  an  improved  understanding  of  the  effects 
of  free-stream  turbulence  in  a  wide  range  of  engineering  problems.  In  a 
modern  gas  turbine  engine,  for  instance,  the  gas  flow  exiting  the  combustor 
contain  high  levels  of  turbulence  which  cause  significant  enhancement  of  heat 
transfer  to  the  downstream  turbine  blades  (Goldstein,  2001),  the  effect  being 
most  severe  in  the  stagnation  point  region  near  blade  leading  edge.  Efforts 
to  improve  the  thermal  eflSciency  and  reliability  of  the  blade  cooling  system 
hinges  critically  upon  an  accurate  prediction  of  heat  transfer  in  the  presence 
of  free-stream  turbulence.  Stagnation  point  flow  also  plays  important  roles 
in  other  industrial  applications  such  as  material  processing  and  electronics 
cooling  (Nakayama,  1995). 

Over  the  years,  a  number  of  experiments  have  investigated  the  heat  trans¬ 
fer  enhancement  over  its  laminar  value  in  stagnation  point  flows  in  the  pres¬ 
ence  of  free-stream  turbulence,  see  Kestin  (1966);  Sadeh  and  Brauer  (1980); 
Van  Fossen,  Simoneau,  and  Ching  (1995);  and  Ames  (1997).  The  turbulence 
intensity,  length  scale  and  the  mean  flow  Reynolds  number  were  shown  to 
be  the  most  important  parameters  in  determining  the  turbulent  heat  trans¬ 
fer  rate.  Typically,  the  heat  transfer  enhancement  was  found  to  increase 
with  increased  Reynolds  number  and  turbulence  intensity,  but  decrease  with 
increased  turbulence  length  scale.  Semi-empirical  correlations  have  been  pro¬ 
posed  to  predict  the  heat  transfer  enhancement,  see  for  example  Smith  and 
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Kuethe  (1966);  Van  Fossen  et  al.  (1995);  and  Ames  (1997).  Numerical  sim¬ 
ulations  have  also  been  performed  to  study  the  detailed  interaction  between 
free-str^m  turbulence  and  a  stagnation  point  boundary  layer.  Spalart  (1989) 
ound  that  out  of  initial  white-noise  disturbances  in  a  swept  Hiemenz  bound- 
ary  ayer,  the  most  unstable  disturbance-mode  is  the  one  that  has  the  same 
^milarity  form  as  the  mean  Hiemenz  flow,  i.e.  the  streamwise  velocity  is  a 
linear  function  of  the  streamwise  coordinate  x  ,  and  the  wall  normal  velocity 
IS  independent  of  x.  The  flow  structures  induced  by  free-stream  turbulence  in 
a  stagnation  repon  axe  found  to  be  qualitatively  similar  to  those  induced  by 
upstream  organized  disturbances  (Xiong  and  Lele,  2001).  The  importance  of 

ir^l  f  has  been  shown  recently  through  direct  numer- 

(2W3)  "  ^  turbulent  stagnation  point  flow  by  Bae,  Lele,  and  Sung 

Theoretically,  the  effect  of  temporal  modulation  of  free-stream  velocity 
was  first  studied  by  Lighthill  (1954)  who  obtained  the  Stokes-layer  correc¬ 
tions  for  skin  friction  and  heat  transfer  for  a  two-dimensional  pulsating  mean 
flow  about  a  cylinder^  The  steady  streaming  (second  order  alteration^o  the 
ean  flow  owing  to  the  Reynolds  stresses)  in  an  oscillatory  Hiemenz  bound¬ 
ary  I^erw^  further  examined  by  Grosch  and  Salwen  (1982)  and  Merchant 

the"^  emphasis  is  on  finding  similarity  solutions  and 

"  to  skin  friction.  The  enhancement  of  heat  transfer  in  a 
perturbed  Hiemenz  boundary  layer  was  also  conceived  as  a  consequence  of 

wZlltnU  ‘he  incoming  disturbances.  However,  Kestin  and 

Wood  (1970)  found,  and  later  much  clarified  by  Wilson  and  Gladwell  (1978) 
hat  the  two-dimensional  Hiemenz  boundary  layer  is  always  linearly  stable 

f  disturbances.  The  nonlinear  instability 

r  'f.^  hyj.yell  and  Huerre  (1985)  who  showed  that  if  the  level  of  the 

^^rtain  threshold  value,  Hiemenz  flow  can  be 
ll^e  bo  n^’  .!  '"®^^hility  for  the  more  general  attachment- 

ind  Thl^r"  investigated  by  Lin  and  Malik  (1996) 

eued  ^^‘^hovin  (1979),  in  a  comprehensive  review,  ar- 

to^P  f  enhancement  of  heat  transfer  is  more  likely  a  forced  response 
itv  4p*^  disturbance  rather  than  a  result  of  internal  flow  instabil- 
ity^The  flow  visualizations  by  Nagib  and  Hodson  (1978)  and  Botcher  and 

vTjr'V'  TV*™"*'/  This  view  poi„rw^  ad 

ocated  earlier  by  Sutera  (1965)  who  analyzed  the  amplification  effect  of  the 

memi  flow  on  the  incoming  disturbances  and  linked  them  to  the  sensitivity  of 

heat  transfer  to  upstream  vortical  disturbances.  By  generalizing  the  classical 
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rapid  distortion  theory  (RDT)  (Batchelor  and  Proudman,  1954),  Hunt  (1973) 
analyzed  the  second  order  moment  of  the  turbulent  velocity  field  when  the 
free-stream  turbulence,  of  either  very  large  or  very  small  integral  scales,  im¬ 
pinges  onto  a  circular  cylinder,  but  the  heat  transfer  between  the  fluid  and  the 
cylinder  was  not  considered.  Dhanak  and  Stuart  (1995)  showed  that,  in  the 
forward  stagnation  region  of  any  two-dimensional  body,  the  viscous  bound¬ 
ary  layer  can  support  a  substructure  of  counter-rotating  streamwise  eddies 
when  there  exists  weak  cross-stream  vorticity  in  the  external  flow.  Kerr  and 
Bold  (1994)  obtained  a  family  of  strained  periodic  vortex  array  embedded  in 
an  inviscid  two-dimensional  stagnation  point  flow.  Andreotti,  Douady,  and 
Couder  (2001)  have  recently  used  these  vortices  as  a  model  to  study  experi¬ 
mentally  the  dynamics  of  interaction  between  strain  and  vorticity.  Although 
much  progress  has  occurred  in  understanding  the  heat  transfer  augmenta¬ 
tion  mechanism,  there  is  however  no,  to  the  authors’s  knowledge,  analytical 
solution  to  the  impinging  vortical  disturbances  in  a  viscous  Hiemenz  bound¬ 
ary  layer.  Theoretical  analysis  on  the  efiects  of  length  scale,  intensity  and 
frequency  of  the  impinging  disturbance  and  the  associated  heat  transfer  has 
been  lacking.  This  is  partly  the  reason  that  the  prediction  of  stagnation  point 
heat  transfer  in  the  presence  of  free-stream  turbulence  has  largely  remained 
empirical.  In  this  paper,  the  distortion  of  unsteady  three-dimensional  dis¬ 
turbances  in  a  two-dimensional  stagnation  point  flow  is  investigated  using 
theoretical  analysis  and  numerical  solutions.  Its  objective  is  to  gain  quanti¬ 
tative  understanding  of  the  heat  transfer  augmentation  mechanism,  particu¬ 
larly  its  dependence  on  disturbance  parameters,  e.g.  length  scale,  intensity, 
and  frequency,  as  a  way  to  improve  the  prediction  of  turbulence  effects  in 
this  technologically  important  flow. 

The  paper  is  organized  as  follows.  The  governing  equations  for  the  mean 
flow  and  the  disturbances  are  formulated  in  §2,  followed  by  a  discussion  of  the 
length  and  velocity  scales  associated  with  the  disturbances.  The  numerical 
solutions  of  the  nonlinear  disturbance  equations,  showing  the  characteristics 
of  the  disturbance  development,  are  presented  in  §3.  Analysis  based  on  a 
linear  vortex  dynamics  is  pursued  in  §4  to  derive  the  dependence  of  vorticity 
amplification  on  the  disturbance  length  and  time  scales.  In  §5,  the  asymptotic 
behavior  for  large  scale  and  low  frequency  disturbance  is  discussed,  along 
with  its  implications  for  the  wall  heat  transfer.  By  superposing  different 
modes  of  upstream  disturbance,  the  analysis  is  extended  in  §6  to  treat  the 
case  of  homogeneous  isotropic  free-stream  turbulence.  A  new  heat  transfer 
scaling  correlation  based  on  the  turbulence  intensity,  integral  length  scale 
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and  Reynolds  number  is  proposed  and  compared  with  recent  experimental 
measurements.  The  conclusions  of  the  present  study  and  a  discussion  is  given 
in  §7. 


6.2  Governing  equations 


We  consider  unsteady,  incompressible,  viscous  flow  with  constant  fluid  prop¬ 
erties  in  the  forward  stagnation  region  of  an  arbitrary  two-dimensional  bluff 
body  shown  in  figure  6.1a.  The  coordinate  axes  are  chosen  as  follows:  x  is 
parallel  to  the  body  and  normal  to  the  attachment  line,  y  along  the  free- 
stream  away  from  the  body,  and  z  along  the  attachment  line.  The  mean  flow 
around  the  bluff  body  is  assumed  to  be  steady  and  two-dimensional  but  the 
incoming  disturbances  are  three-dimensional  and  may  vary  with  time.  The 
length  scale  of  the  disturbances  is  assumed  to  be  large  compared  with  the 
boundary  layer  thickness  but  much  smaller  than  the  diameter  of  curvature  at 
the  stagnation  point,  hence  the  mean  flow  in  this  region  is  well  modelled  by  a 
plane  Hiemenz  boundary  layer  flow,  see  for  example  in  Batchelor  (1967).  In¬ 
deed,  Wilson  and  Gladwell  (1978)  showed  that  as  Reynolds  number  Re  — >  oo, 
the  laminar  flow  in  the  stagnation  region  of  any  two-dimensional  bluff  body 
is  reduced  to  a  plane  Hiemenz  flow  problem.  Exploiting  this  reduction,  the 
present  analysis  will  be  focused  on  the  Hiemenz  boundary  layer  flow  in  the 
presence  of  upstream  disturbances  as  shown  in  figure  6.1b.  By  relating  the 
strain  rate  in  Hiemenz  solution  to  the  free-stream  velocity  and  the  diameter 
of  curvature  at  the  stagnation  point,  the  present  analysis  applies  to  a  general 
two-dimensional  bluff  body. 

When  the  characteristic  length  scale  Iq  and  velocity  scale  Uo  of  the  Hiemenz 
boundary  layer,  deflned  as 


h  =  \/u*/A*, 


Uq  =  y/i/*A* 


(6.1) 


are  used  to  nondimensionalize  the  coordinates  and  the  flow  variables,  we 
have 


r),  C) 

(u,  V,  w) 


(L.  El  £.) 

Wo’  /o’  lo^ 

(—  E1 
\  >  >  ) 
Uo  Uo  Uo 


(6.2a) 

(6.2b) 


P  = 


(6.2c) 
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(6.2d) 


where  the  superscript  *  is  used  hereafter  to  denote  the  dimensional  quantities, 
V*  is  the  kinematic  viscosity  and  vl*  is  the  strain  rate  of  the  external  potential 
flow  in  the  Hiemenz  solution,  p*  and  p*  are  the  pressure  and  density.  The 
wall  is  assumed  to  be  isothermal  with  temperature  and  the  upstream  flow 
temperature  is  T^. 

The  flow  field  {a,p,6}  is  assumed  to  consist  of  a  plane  stagnation  point 
flow  {U,P,  0}  and  a  unsteady  disturbance  field  {u,p,  0},  i.e. 

u  =  u+u,  p  =  p+p,  e  =  e  +  e  (e.s) 

Following  Batchelor  (1967),  the  mean  velocity  U  of  Hiemenz  flow  may  be 
expressed  as 

V  =  {U,V,W)  =  (<l/^,-(f>,0).  (6.4) 

Together  with  the  mean  temperature  0,  they  satisfy  the  following  Hiemenz 
equations 

4>"'  +  <1><I>"  +  l-(f>'^  =  0  (6.5a) 

0"  +  Fr<^0'  =  O.  (6.5b) 

where  ^  is  only  a  function  of  rj  and '  denotes  d/dr].  Pr  is  the  Prandtl  number. 
The  boundary  conditions  for  (f)  and  0  are  given  by 

<^(0)  =  .^'(0)  =  ^'(oo)  =  1  (6.6a) 

0(0)  =  0,  0(oo)  =  1  (6.6b) 

The  general  governing  equations  for  the  perturbation  field,  following  from 
(6.1)  and  (6.3),  can  be  written  as 

V  •  u  =  0  (6.7a) 

5iu  +  u-Vu  +  U-Vu  +  u-VU  =  -Vp  +  V2u  (6.7b) 

5*0  -I-  u  •  +  U  •  V0  +  u  •  V0  =  (6.7c) 

Pr 

In  this  paper,  perturbations  of  the  form 

u  =  viViCit),  w{i],C^t)},  P  =  p{v,C,t),  0  =  9{r],C,t) 

(6.8) 


148 


are  considered;  the  perturbation  variables  have  no  ^  dependence  except  the 
e-component  of  u.  As  noted  by  Spalart  (1989),  this  is  a  good  approximation 
to  the  flow  field  in  the  neighborhood  of  the  stagnation  ^  plane.  Steady 
perturbations  of  this  form  were  also  used  by  Sutera  (1965).  Substituting 
(6.8)  into  (6.7)  yields  the  governing  equations  for  the  perturbation  field: 


u  +  w'  +  d^w  =  0 


(6.9a) 


dtu  +  (u^  +  vti'  +  wd^u)  -(fn/  +  v<f/'  +  2u<f/  =  u"  +  d^u 
dtv  +  {vv'  +  wd^v)  -  {(fwy  =  -p'  +  (u"  +  d^v) 
dtw  +  (W  +  wd^w)  -  <fnv'  =  -d^p  +  {w"  +  d^w) 


(6.9b) 

(6.9c) 

(6.9d) 


dt9  +  (v^  +  wd^9)  +  v&-(f>9'=:  ^(^'  +  ^0) 


(6.9e) 


Note  here  the  prime  '  denotes  the  partial  derivative  The  disturbance 
vorticity  may  be  conveniently  found  as 


u}^  =  w'-diV,  or,  u}i  =  -^u'  (6.10) 

As  will  be  seen  in  subsequent  discussion,  the  disturbance  vorticity  in  ^  di¬ 
rection  Wf,  which  is  subject  to  stretching  by  the  mean  diverging  flow,  plays 
a  central  role  in  describing  the  disturbance  evolution  in  this  type  of  flow. 
Hence  the  governing  equation  for  streamwise  vorticity,  hereafter  denoted  by 
u,  is  written  as 


dtu  -  (a;"  -f  ^o;)  -  {tfxvY  =  -{yuf  -  d^{wu})  (6.11) 

To  seek  solutions  which  are  periodic  in  time  t  (or  steady)  and  periodic 

in  the  spanwise  direction  C,  we  expand  u,v,w,o)  and  ^  in  a  double  Fourier 
series: 


00 

=  Aj,  ^  UmniTi)exp{i{Tnaot  +  nkoO} +v'q{tj)  (6.12a) 

m,n=l 

oo 

HvXyt)  =  Ap  Vmn(v)exp{i{Tnaot  +  nko(:)} -vo{Tj)  (6.12b) 

m,n=l 

oo 

=  Ap  {nko)~^Wmn{v)exp{i{maot  +  nkoQ}  (6.12c) 

m,n=l 
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oo 

^mn{v)  exp{z(m(7ot  +  nfcoC)}  {6.12d) 

m,n=l 
oo 

=  Ap'Yl  ^mn(^)exp{i(7n<7oi  +  nA:oC)}+^o{»?)  (6.12e) 

m,n=l 

where  Ap  is  the  perturbation  amplitude;  fco  =  k*lQ  is  the  fundamental  wavenum¬ 
ber  in  spanwise  direction  and  oq  =  a* /A*  the  fundamental  frequency.  The 
functions  Vo{rf)  and  0o(?7)  represent  the  nonzero  spanwise  averages  of  the 
perturbation  velocity  and  temperature,  i.e.  the  modification  to  the  mean 
flow  profiles  due  to  the  nonlinear  interaction  among  the  disturbance  modes. 
Additionally,  the  normal  derivative  of  ^o(^)  at  the  wall  gives  the  spanwise 
averaged  heat  transfer  enhancement. 

Substituting  the  expansions  (6.12)  into  the  continuity  equation  (6.9a) 
yields 

«mn  +  +  *  '^mn  =  0,  (6. 13) 

and  from  the  definition  of  Wmn?  it  follows  that 

i  n^kl  Vmn  -  Wmn  +  ^mn  =  0.  (6-14) 


Similarly,  the  governing  equations  for  streamwise  velocity  Umn,  vorticity  oJmn 
and  temperature  dmn  reduce  to 


“mn  -  (2^*  +  ^^^0  +  i'm<^o)Umn  +  <{>  «mn  “  ^mn  =  M){u,  V,  w)  (6.15) 

^mn  -  +  imao)u}mn  +  {(f> t^mnY  =Afl{v,W,u)  (6.16) 

-  {ri^kl  +  imPrao)  ^mn  +  -  Prvmn0  =  A/’2(u,  w,  9)  (6.17) 

The  equations  for  Uo(^)  and  0o(^)  can  also  be  derived  as 

Vq  +  <fyvQ  -  2<f/vQ  +  <f/'vo  =  A/3(«,  v)  (6.18) 

Pr(f>  9'o  =  M,(v,  w,  9)  (6.19) 

In  the  above  equations,  A/i’s  are  the  nonlinear  terms  given  by  the  following 
expression: 


A/i)  = 


flE 

2 


{^pg'^m-p,n-g  +  UpgUm+p,n+g  +  ^P9^m-p,n-g 

p,g=l 


+  VpqU^^p.^^g 
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+  (tI  +  9)'^i^j>^Wf7i+p,n+g]  } 

(6.20a) 

A/i 

A  ^ 

2  ^  ^  “1" 
p,g=l 

m. 

+  ^  ]}^pq^Tn-p,n-q  +  '^pq^m-^p,n-\-q]  } 

(6.20b) 

A4 

Aj,Pr^ 

~  2  +  ^PqKn+p,n+q 

P.?=l 

+“[(**  ~  Q)‘^pq^tn-p,n-q  +  («  +  Q)'^p^mA^p,n-A-^  } 

(6.20c) 

A/3 

00  1 

p.9=l 

(6.20d) 

Na 

A^Pr  °° 

■“  2  ^  ^  ^^{^M^pq  "I"  ^“^pq^pq^ 

P,9=l 

(6.20e) 

where"  stands  for  the  complex  conjugate.  Suppose  that  unsteady  disturbance 
vorticity  is  introduced  upstream  at  ?/  =  /fo  »  1  by  superimposing  a  simple 
sinusoidal  variation  with  the  amplitude  Ap,  at  the  fundamental  spanwise 
wavenumber  ko,  and  frequency  (Tq  on  the  mean  velocity  The  disturbance 
boundary  conditions  a.t  t)  =  Hq  are 

*'11  =  1,  Vmn  =  0  for  m,n^  1,  Umn  =  Wmn  =  ^mn  =  0  (6.21) 

At  the  wall  17  =  0,  no  slip  and  isothermal  boundary  conditions  are  enforced 
for  the  velocities  and  temperature,  thus 


^mn  —  Wmn  —  Wmn  —  ^mn  —  0  (6.22) 

Once  the  disturbance  velocity  v^n  ,  Wmn  and  temperature  9mn  are  obtained, 
the  relative  heat  transfer  enhancement  A/i  over  its  undisturbed  mean  value 
h  can  be  found  by  solving  (6.19) 


h  e'(O) 


(6.23) 


It  is  convenient  to  introduce  new  length  and  velocity  scales  besides  those 
in  (6.1).  By  the  assumed  spanwise  periodicity,  a  natural  choice  for  the  dis¬ 
turbance  length  and  velocity  scales  are 


Id  =  lA*,  Ud  =  vk* 


(6.24) 
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It  is  observed  below  that  the  square  of  the  ratio  between  the  disturbance 
length  scale  la  and  the  Hiemenz  boundary  layer  scale  lo,  represented  by  a 
dimensionless  parameter  A  see  Kerr  and  Dold  (1994), 

A  =  (—\^  —  ^  —  _L 

kl 

is  a  critical  parameter  in  determining  the  evolution  of  the  disturbances. 
When  A  is  large,  the  distortion  of  the  upstream  disturbances  is  mainly  due 
to  the  mean  flow  straining  effect;  the  nonlinear  interaction  among  the  dis¬ 
turbance  modes  is  of  higher  order.  This  is  similar  to  the  cases  treated  by 
the  traditional  linear  RDT.  But  in  the  presence  of  viscosity,  the  nonslip  wall 
introduces  the  viscous  effect  in  its  vicinity  and  must  be  included  in  the  for¬ 
mulation  both  for  the  mean  flow  and  the  disturbances.  Interestingly,  A  can 
also  be  interpreted  as  the  time  scale  ratio  between  the  disturbance  turn-over 
time  Idfud  and  the  mean  flow  straining  time  \/A*,  i.e. 

A  =  {^)A*  (6.26) 

Ud 

For  different  values  of  A  as  well  as  Ap,  the  numerical  solutions  to  the  flow 
problem  posed  here  are  presented  in  the  next  section. 

6.3  Numerical  Results 

The  system  of  equations  (6.13)  -  (6.19),  with  boundary  conditions  (6.21)- 
(6.22),  forms  a  second  order  boundary- value  problem  driven  by  an  inhomo¬ 
geneous  boundary  condition  away  from  the  wall.  Numerically,  they  can  be 
readily  solved  using  the  over-relaxation  method,  described  for  example  in 
Isaacson  and  Keller  (1993).  Fourth  order  finite  difference  scheme  is  adopted 
to  approximate  the  spatial  derivatives  and  different  numbers  of  grid  points 
are  used  to  obtain  the  grid  independent  solutions.  For  the  nonlinear  calcula¬ 
tions,  a  total  number  of  modes  resulting  from  a  truncation  at  m  =  6,  n  =  6 
of  the  double  Fourier  series  in  (6.12)  are  found  sufficient  for  the  solutions  to 
converge.  The  details  of  the  numerical  method  and  convergence  study  can 
be  found  in  Xiong  (2004). 

The  case  of  steady  disturbance,  i.e.  <7o  =  0,  is  discussed  first;  this  will 
help  clarify  the  dependence  of  the  disturbance  evolution  on  the  length  scale 
ratio  A.  In  all  computations,  the  Prandtl  number  is  taken  as  Pr  =  0.71, 
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and  the  disturbance  amplitude  Ap  is  ten  percent  of  the  mean  flow  velocity 
0  at  the  inflow  boundary  1/  =  Hq,  i.e.  Ap  =  0.1<f>{Ho)  and  Ho  =  18.  Figure 
6.2  shows  the  profiles  of  the  base  state  velocity  and  temperature  (without 
any  free-stream  disturbance).  The  velocity  and  temperature  boundary  layer 
thicknesses  (defined  as  the  location  at  which  99%  of  the  external  value  is 
reached)  are  around  rf  =  2.4  and  t/  =  4,  respectively.  Figure  6.3  shows  the 
streamline  patterns  of  the  perturbed  flow  in  the  stagnation  plane  for  different 
values  of  A.  The  horizontal  axis  spans  one  spanwise  wavelength  for  each  case. 

As  can  be  seen,  initially  at  A  =  1.1,  the  perturbed  streamlines  converge 
toward  a  free  stagnation  point  at  the  symmetry  plane  z  -0.  At  A  /K/  0(2), 
a  pair  of  counter  rotating  vortices  start  to  form  at  the  edge  of  the  boundary 
layer.  The  strength  of  these  vortices  increases  rapidly  with  the  increasing  A 
and  attains  its  maximum  at  A  =  4.  When  A  increases  further,  the  vortex 
strength  slowly  decreases.  At  A  =  32  the  vortices  disappear  and  are  replaced 
again  by  a  sink  type  free  stagnation  point.  For  even  larger  A  =  64,  the 
mean  flow  dominates,  and  the  free  stagnation  point  also  disappears.  Qual¬ 
itatively,  these  streamline  patterns  can  be  first  classified  into  two  groups 
depending  whether  a  free  stagnation  point  (FSP)  is  present.  The  stagna¬ 
tion  point  emerges  when  the  disturbance  wall  normal  velocity  v  exceeds  the 
mean  flow  Since  at  the  wall,  <f)  =  <f>'  =  v  =  =  0,  this  can  occur  when 

k  (0)1  >  \<t>  (0)|.  Furthermore,  depending  on  the  direction  of  the  spanwise 
velocity  in  the  neighborhood  of  FSP,  the  resulting  FSP  can  be  either  a  sink 
point  when  the  spanwise  velocity  w  points  inward,  or  a  saddle  point  — 
when  it  points  outward.  Since  by  symmetry  tu  =  0  at  z  =  0,  the  direction  of 
w  is  determined  by  From  continuity  equation,  this  is  in  turn  determined 
which  the  vertical  velocity  tends  to  zero  at  the  free  stag¬ 
nation  point.  When  ^  <  0,  it  becomes  a  saddle  point  and  the  streamlines 
emanating  from  the  stagnation  point  eventually  form  the  spiral  vortices.  For 
the  present  case  with  fixed  Ap,  the  free  stagnation  point  emerges  for  A  >  1, 
and  transition  from  a  sink  point  to  a  saddle  point  type  occurs  at  A  =  1.6.  At 
A  =  30,  the  saddle  type  stagnation  point  changes  back  to  a  sink  point  type, 
causing  the  vortices  to  disappear.  At  A  >  60,  the  free  stagnation  point  also 
disappears,  and  the  perturbed  flow  becomes  unidirectional  in  wall  normal 
direction. 

In  addition  to  the  change  of  A,  the  effect  of  the  disturbance  amplitude 
Ap  on  the  vortex  formation  is  shown  in  figure  6.4  where  Ap  varies  from  2 
to  15  percent  of  ^(Hq)  for  fixed  A  =  4.  In  this  case,  the  free  stagnation 
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Amplitude  Ap  f  <l>{Ho) 

Velcoity 

Vorticity  |a;tu| 

Heat  transfer  Ah  /  h 

0.02 

0.9321 

1.4953 

0.0100 

0.05 

2.5558 

3.7457 

0.0620 

0.07 

3.8928 

5.2490 

0.1206 

0.10 

6.4999 

7.5256 

0.2406 

0.12 

8.8094 

9.0672 

0.3368 

0.15 

13.729 

11.410 

0.4918 

Table  6.1:  The  maximum  wall  normal  disturbance  velocity  |u„,oa:l,  vorticity 
at  the  wall  and  the  heat  transfer  enhancement  Ah/h  for  different  Ap  at 
A  =  4,Fo  =  18. 

point  first  emerges  at  Ap  =  0.03(^(i/o)>  the  transition  from  the  sink 
type  to  the  saddle  type  takes  place  at  Ap  =  0.05^(//o)-  The  latter  can  be 
considered  as  a  threshold  for  Ap  since  only  beyond  this  value  the  counter¬ 
rotating  vortices  become  possible  in  the  streamline  pattern.  Moreover,  figure 
6.4  also  reveals  that  although  the  strength  of  the  vortices  increases  with  the 
increasing  Ap,  the  overall  flow  pattern  remains  qualitatively  similar  once  the 
threshold  is  reached.  The  existence  of  a  threshold  in  disturbance  amplitude 
is  also  consistent  with  earlier  observations  by  Nagib  and  Hodson  (1978)  on 
the  formation  of  vortex  pair  at  the  stagnation  region  of  a  bluff  body  subject 
to  the  impingement  of  wakes.  Similar  flow  patterns  as  those  shown  in  figure 
6.3  and  6.4  also  emerge  in  aforementioned  experiments  (Nagib  and  Hodson, 
1978;  Botcher  and  Wedemeyer,  1989)  as  well  as  numerical  simulations  (Xiong 
and  Lele,  2001;  Bae  et  al.,  2003).  The  quantitative  characterization  of  the 
flow  fields  in  figure  6.4  in  terms  of  the  maximum  wall  normal  disturbance 
velocity  jumoxlj  wall  vorticity  |a;,„|  and  heat  transfer  enhancement  Ah/h  are 
summarized  in  Table  6.1. 

Figure  6.5  shows  the  profiles  of  the  fundamental  mode  and  higher  har¬ 
monics  at  A  =  4  for  the  disturbance  velocity,  vorticity  and  temperature.  The 
disturbance  v  velocity  and  vorticity  are  found  to  be  amplified  before  reach¬ 
ing  the  boundary  layer.  Compared  to  the  fundamental  mode  m  =  l,n  =  1, 
the  amplitudes  of  the  higher  harmonics  are  typically  sjnall,  except  for  tem¬ 
perature  where  the  mean  temperature  modification  Oq  attains  an  amplitude 
similar  to  ^n.  Of  the  three  components  of  the  disturbance  velocity,  Umn  is 
found  to  be  typically  one  order  of  magnitude  smaller  than  Vmn  and  Wmn-  The 
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corresponding  contours  for  the  vorticity  and  temperature  are  plotted  in  fig¬ 
ure  6.6  and  figure  6.7.  The  incoming  vorticity  is  amplified  by  a  factor  greater 
than  5  due  to  stretching  and  thus  large  amplitude  vorticity  with  opposite 
sign  is  induced  within  a  thin  near-wall  region  to  satisfy  the  no  slip  boundary 
condition.  The  temperature  contours  are  also  significantly  modified  by  the 
disturbance  velocity.  The  upward  velocity  causes  the  local  thermal  boundary 
layer  to  become  thicker,  while  the  downward  velocity  makes  it  thinner.  The 
net  effect  on  the  spanwise  averaged  heat  transfer  depends  on  the  strength 
and  distribution  of  these  thickened  and  thinned  regions. 

In  figure  6.8,  the  relative  heat  transfer  enhancement  is  plotted 

as  a  function  of  A.  For  small  A,  or  small  disturbance  wavelength,  Ah/h 
increases  rapidly  as  A  increases.  But  for  large  A  ,  it  decrease  slowly  with 
the  increase  of  A.  A  peak  value  is  found  around  A  =  4,  indicating  that  an 
optimum  disturbance  length  scale  exists  at  about  five  times  the  boundary 
layer  thickness  and  produces  maximum  heat  transfer  enhancement.  Also 
shown  in  figure  6.8  is  the  maximum  value  of  the  fundamental  mode  of  the 
disturbance  wall  normal  velocity  vn-  It  shows  a  similar  trend  with  A  and 
optimum  amplification  at  nearly  the  same  value  of  A.  However,  for  large  A, 
Vmax  ~  A“^/*  while  Ah/h  ~  A~^  Finally,  in  figure  6.9  the  amplitude  of  the 
vorticity  at  the  wall  is  plotted  as  a  function  of  A.  Besides  a  similar  peak 
location  around  A  =  4,  the  striking  feature  of  the  wall  vorticity  is  that  it 
approaches  a  constant  which  depends  only  on  Hq  and  Ap  as  A  becomes  large. 

The  dependence  of  flow  characteristics  on  the  length  scale  is  due  to  the 
competition  between  the  vortex  stretching  and  the  viscous  dissipation.  The 
convective  heat  transfer  or  more  generally  passive  scalar  transport  is  a  direct 
consequence  of  the  amplified  velocity  disturbances,  and  may  in  fact  be  re¬ 
garded  as  an  indication  of  how  significantly  the  flow  near  the  wall  has  been 
modified.  Based  on  this  point  of  view,  the  vorticity  equation  (6.16),  compris- 
ing  the  mean  flow  stretching,  viscous  dissipation  and  nonlinear  interaction 
effects,  shall  be  the  starting  point  for  analyzing  the  disturbance  evolution. 
From  the  profiles  of  in  figure  6.5,  the  nonlinear  interaction  is  shown  to 
be  relatively  weak,  hence  in  the  subsequent  analysis  the  nonlinear  terms  in 
(6.16)  will  be  neglected.  A  similar  rationale  is  also  the  basis  for  the  linear 
RDT  which  has  been  successfully  applied  to  this  type  of  flow.  However,  un¬ 
like  the  purely  inviscid  interaction  considered  in  R,DT,  in  the  present  problem 
viscosity  exerts  a  significant  influence  upon  the  disturbance  throughout  the 
whole  domain,  even  when  the  mean  flow  can  be  treated  largely  as  inviscid  in 
the  region  outside  the  boundary  layer.  Indeed,  it  is  the  balance  between  the 
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vortex  stretching  and  viscous  diffusion  that  produces  the  maximum  vorticity 
near  the  edge  of  boundary  layer.  Therefore,  analysis  of  the  vorticity  dy¬ 
namics  helps  to  understand  the  basic  mechanism  governing  the  flow  and  the 
associated  scalar  transport.  This  is  carried  out  in  detail  in  the  subsequent 
sections. 


6.4  Linear  Vortex  Dynamics 

The  primary  goal  of  the  analysis  in  this  section  is  to  find  an  asymptotic 
expression  describing  the  evolution  of  the  large  scale,  low  frequency  distur¬ 
bances  in  the  Hiemenz  boundary  layer.  The  solution  is  sought  by  first  ex¬ 
pressing  the  vorticity  in  a  series  expansion  based  on  the  large  length  scale 
(A  »  1)  and  low  frequency((7o  <C  1).  This  solution  is  formally  valid  in  the 
entire  spatial  domain,  but  insufficient  to  describe  the  vorticity  evolution  in 
an  explicitly  way  owing  to  the  lack  of  closed  form  expression  for  <j>.  On  the 
other  hand,  by  exploiting  the  particularly  simple  form  taken  by  0  in  the 
region  outside  the  Hiemenz  boundary  layer,  a  closed  form  solution  for  the 
vorticity  can  be  found  for  any  arbitrary  A  and  ao  in  that  region.  These  two 
solutions  are  required  to  match  in  the  region  outside  the  Hiemenz  boundary 
layer  where  they  are  both  valid.  Thus,  an  explicit  composite  asymptotic  so¬ 
lution  can  be  formulated  which  describes  the  evolution  of  the  large  scale,  low 
frequency  vortical  disturbances,  and  forms  the  basis  of  further  analysis. 

6.4.1  Series  expansion 

As  shown  by  the  numerical  results  in  figure  6.5,  the  disturbance  typically 
reaches  its  maximum  amplitude  at  the  edge  of  the  boundary  layer  before 
it  decays  .  By  (6.16)  and  (6.25),  the  linearized  governing  equation  for  comn 
takes  the  form 


2 

”  (y  +  *  =  0  (6.27) 

Here  only  the  case  of  large  A  1  and  small  ctq  <C  1  will  be  considered.  These 
limits  correspond  to  the  situation  where  the  upstream  disturbance  is  of  large 
scale  and  low  frequency  relative  to  the  Hiemenz  boundary  layer  scales.  This 
regime  is  prototypical  for  the  free  stream  turbulence  impinging  on  gas  turbine 
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blades.  A  small  parameter  c  may  be  defined  as: 


n‘  .  ,  , 

^  =  y  +  I  moQ,  |c|  <  1  (e.28) 

and  a  series  solution  of  is  sought  in  the  powers  of  e  as  a  regular  pertur¬ 
bation: 

u;mn  =  a;®„-t-€a;^„  +  ...  (6.29) 

The  equation  at  zeroth  order  becomes 

^mnY  =  0  (6.30) 

whose  general  solution  can  be  found  as 


<^L{v)  =  e-*  +  F® e-*  dr)'. 

Jo 

At  first  order  the  equation  is 

«„r  +  ((»<„)'  =  <n 

and  the  solution  can  be  expressed  using  as 


(6.31) 


(6.32) 


<»(>?)=  E'e-*  +  F'e-» 

•'0  Jo  Jo 

where  $  is  defined  as 

and  F®,F\  F®,Fi  are  arbitrary  constants.  The  higher  order  terms  can  be 
computed  similarly.  The  series  expansion  (6.29)  is  valid  in  the  entire  spatial 
domain  from  the  wall  to  the  inflow  but,  as  such,  of  limited  use  because  it 
involves  unknown  constants.  In  next  section,  the  exact  solution  of  (6.27)  in 
a  region  outside  the  Hiemenz  boundary  layer  is  obtained  for  arbitrary  A  and 
<7o-  By  examining  the  characteristics  of  the  exact  solution  for  the  case  of 
A  >  1  and  (To  <C  1,  the  unknown  constants  in  (6.31)  can  be  determined. 


f  <f>{v')dT}' 

Jo 


(6.33) 

(6.34) 
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6.4.2  Outside  the  boundary  layer 

Outside  the  Hiemenz  boundary  layer,  the  relevant  length  and  velocity  scales 
are  determined  by  the  incoming  disturbance.  To  facilitate  the  subsequent 
analysis,  the  vorticity  equation  (6.27)  is  first  rescaled  by  1^  and  Uj  defined  in 
(6.24) 

w"  „  +  X<f)  +  {X4>  -r?  -i  mXao)  comn  =  0.  (6.35) 

Note  here  both  the  dependent  and  independent  variables  are  rescaled,  and 
the  derivatives  are  now  with  respect  to  the  new  independent  variable  s  = 
y*/ld  =  T)/VX-  In  particular,  the  mean  velocity  profile  (f>{s)  in  (6.35)  is  equal 
to  the  multiplied  by  the  factor  1/y/X  ,  i.e.  <^(s)  =  The 

boundary  conditions  for  the  mean  and  disturbance  flow  remain  unchanged 
under  this  rescaling. 

As  noted  in  figure  6.2,  the  mean  velocity  profile  0  takes  a  simple  irrota- 
tional  form  outside  the  Hiemenz  boundary  layer 

^  s  —  Sd,  (f/  <^1.  (6.36) 

where  Sd  represents  the  displacement  thickness  Sd  of  the  Hiemenz  bound¬ 
ary  layer,  i.e.  Sd  =  Sd/Vx.  In  this  region,  the  main  effect  of  (f>  on  the 
disturbances,  besides  the  convection,  is  to  stretch  the  ^—component  of  the 
incoming  vorticity.  When  the  disturbance  scale  is  relatively  large,  viscosity 
plays  a  less  important  role  and  the  straining  effect  leads  to  an  increase  of  the 
vorticity. 

Denoting  Umn  outside  the  boundary  layer  by  using  (6.36),  (6.35) 

becomes 

+  A  (s  -  Sd)  «„)'  -I-  (A  -  -  2  mAao)  <„  =  0  (6.37) 

On  introducing  a  new  independent  variable  r 

r  =  -^A(s  -  Sdf  (6.38) 

the  vorticity  equation  is  further  transformed  into 

tKJ"  +  (|  -  r){AJ  -  ~  =  0  (6.39) 

where  '  stands  for  d/dr.  This  is  the  confluent  hypergeometric  equation  of 
the  general  form 

xy"  +  (c  —  x)y'  —  ay  =  0  (6.40) 
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whose  solution  may  be  expressed  as 

y  =  Cl  A/ (a;  c;  x)  +  C2  C(c  —  a;  c;  — x)  e® 

for  c  ^  0,  ±1,  ±2  •  - M {a;  c;  x)  and  C(a;  c;  x)  are  the  first  and  second  kind 
of  confluent  hypergeometric  functions  (Abramowitz  and  Stegun,  1970).  In 
the  present  case,  they  correspond  to 


_ X-v?-  imXao 

2’  2A 


y  ^ 

An  integral  representation  for  M(a;c;x)  when  Rec>  Rea>0  is 


(6.41) 


=  rra)rfcL)  -  *r-'^dt 

and  f/(a;c;x)  can  be  expressed  by  M(a;c;x)  through  the  following  expres- 
Sion. 


U(a;  c;x)  = 


sin(c7r)'r(l  +  a-  c)r(c) 


—  X 


i-cAf(l+a-c;  2  -  6;  a:) 


r(a)r(2  -  c) 


} 

(6.42) 


Honco  the  general  solution  for  the  vorticity  becomes 

1;  L,  (6.43) 

For  A  <  1,  the  solution  simply  decays  from  its  upstream  value,  and  for  A  =  1 
it  becomes  a  constant.  So  in  what  follows  only  the  case  of  A  >  1  is  considered. 
First  writing  s  in  terms  of  the  original  variable  t},  and  by  the  definition  of  r 
it  follows  ’ 

n  =  y/Xs,  ^  =  -— =-|-  (6.44) 

where  the  shifted  coordinate  Tj  is  defined  as 

f}  =  r)-5d.  (6.45) 

After  some  simplification,  the  general  solution  for  the  vorticity  outside  the 
boundary  layer  may  be  expressed  as 


Af  ( 


2A 


2A 


~*“2^’  2’ 


(6.46) 


where  Cmn,  ^mn  are  arbitrary  constants.  A  typical  solution  for  A  =  4, 
<ro  =  0  composed  of  two  confluent  hypergeometric  functions  (  with  C, 
Dfnn  specified  by  (6.61)  in  section  6.4.4)  is  illustrated  in  figure  6.10 


'mn? 
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6.4.3  Composite  solution 

To  facilitate  further  analysis,  the  series  solutions  (6.29)-(6.33)  is  combined 
with  the  solution  outside  the  boundary  layer  to  derive  a  composite  asymptotic 
solution.  To  do  so,  the  unknown  coeflBcients  in  the  series  expansion  need  to 
be  determined  first.  Notice  that  as  t)  oo, 

fi^ 

ij>  fj  and  $  ^0  +  -^  (6-47) 

oo  6^ 

where  $o  =  ~  ^  is  an  integral  constant.  This  asymptotic 

representation  is  valid  for  fj  >  fjo,  where  is  a  location  beyond  which 

the  irrotational  flow  applies.  The  zeroth  order  series  solution  (6.31)  may  be 
rewritten  as 

rvo  rn 

+  F^  Uri']  +  F^  >  dr)’.  (6.48) 

Jo  Jtjo 

However  when  rj  — >  oo,  the  first  term  in  (6.48)  vanishes,  so  it  follows 

cOmniv)  ~  f  drf  ~  f  dff.  (6.49) 

Jrto  Jo 

where  in  the  second  expression,  (6.47)  has  been  used  and  the  lower  limit  has 
been  extended  to  zero  using  the  same  decomposition  in  (6.48).  On  the  other 
hand,  the  corresponding  zeroth  order  expansion  of  (6.46)  in  e  becomes 

<„  =  1;  -^)  +  |  -^)v.  (6.50) 

Noting  the  following  identities 

M{a;  a;  z)  =  e*  and  e~^  f  e^dz'  =  zM(l;  — ^)  (6.51) 

Jq  2  2 

and  comparing  (6.49)  and  (6.50)  for  large  tj,  we  obtain 

Dmn  =  (6.52) 

since  the  first  term  in  (6.50)  vanishes  more  rapidly  than  the  second. 

No  appropriate  matching  condition  for  Cmn  can  be  derived  from  the  large 
TJ  asymptotics.  Hence,  it  is  necessary  to  find  the  matching  condition  by 
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considering  the  behavior  of  the  vorticity  close  to  the  wall.  As  shown  in 
figure  6.9,  the  wall  vorticity  plays  a  key  role  in  determining  the  interaction 
between  the  incoming  disturbances  and  the  Hiemenz  boundary  layer.  An 
initial  choice  for  matching  Cmn  S'Od  therefore  seems  to  require  the  wall 
vorticity  obtained  from  (6.46)  and  (6.31)  to  be  equal.  However,  the  vorticity 
at  the  wall  obtained  by  setting  77  =  0  in  (6.46)  is  not  accurate,  since  (6.46) 
is  only  applicable  outside  the  Hiemenz  boundary  layer.  A  more  appropriate 
choice  is  to  require  the  vorticity  obtained  from  the  two  solutions  to  match 
at  the  edge  of  the  Hiemenz  boundary  layer  where  both  solutions  apply.  This 
results  the  following  expression  for 


Cmn  -  [El^  +  e*  dT}]e  **  -  J  dr]  (6.53) 

where  =  ^{S)  -  {S  -  Saf/2.  Further  matching  (6.46)  to  (6.33)  at  higher 
order  requires  evaluation  of  M(a,c,z)  subject  to  the  perturbation  of  a  and 
the  asymptotic  behavior  of  the  last  term  in  (6.33)  for  large  77,  for  which  the 
explicit  expressions  have  not  been  obtained.  Nevertheless,  an  approximate 
solution  of  iOmn,  more  explicit  than  the  series  expansion,  can  still  be  con¬ 
structed  by  a  similar  procedure  to  those  in  matched  asymptotics  (Van  Dyke, 
1975),  i.e.  ’ 

Wmn  ~  +  uC„  -  (6.54) 

here  the  expressions  for  a;«„,  and  are  (6.31),  (6.46)  and  (6.50), 

respectively.  Using  the  matching  condition  (6.52)  and  (6.53),  we  obtain  the 
following  asymptotic  expression  for  the  incoming  vorticity  in  the  case  of 
A  ^  1  and  <to  1  in  the  whole  range  of  77 


Cmn{M{^^-^  -  i^-  1;  -i?.  j 


2A 

/ 

Jo 


2  ’  2 


■)(’!  -  id)  + 


e  *  dr]-^e' 


(6.55) 


Here  the  dependence  of  Umn  as  an  explicit  function  of  A  and  oq,  not  readily 
obtained  by  the  series  expansion  itself,  has  been  retained  in  (6.55),  and  by 
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correcting  the  in  the  near  waJl  region  through  (6.55)  also  extends 
inside  the  Hiemenz  boundary  layer.  In  figure  6.11  and  figure  6.12,  the  com¬ 
parisons  between  the  composite  asymptotic  solution  (6.55)  and  the  numerical 
solution  are  shown  for  A  =  4  and  A  =  36,  respectively.  The  vorticity  pro¬ 
files  have  been  normalized  by  its  value  at  the  inflow  boundary  rj  =  Hq.  The 
numerical  solution  are  obtained  by  solving  both  the  linear  equation  (6.35) 
and  the  fully  nonlinear  equation  (6.16).  The  asymptotic  solution  and  the 
linear  numerical  solution  agree  well  at  A  =  4  and  become  indistinguishable 
at  A  =  36.  Recall  that  the  disturbance  amplitude  is  here  taken  as  10  per¬ 
cent  of  the  mean  firee-stream  velocity,  Ap  =  O.10(/fo)j  and  the  linear  solution 
approaches  the  nonlinear  solution  quite  well  as  A  increases.  This  indicates 
that  the  present  linear  analysis  is  adequate  for  describing  the  characteristics 
of  disturbance  development  in  a  Hiemenz  flow. 

6.4.4  Boundary  conditions  for  the  vorticity 

The  general  composite  solution  in  (6.55)  describes  the  evolution  of  distur¬ 
bance  vorticity  in  the  whole  domain  from  the  inflow  to  the  wall  boundary. 
Now  we  set  to  specify  the  constants  Cmn  and  Dmn  by  the  inflow  and  wall 
boundary  conditions.  In  fact,  it  will  later  become  clear  that  both  Cmn  and 
Dmn  bear  clearly  defined  physical  meanings.  First,  the  initial  disturbance 
vorticity  is  introduced  at  the  inflow  boundary  far  upstream,  i.e. 

Wmn  =  <^mn{H(i)  at  T)  =  Hq  (6.56) 

Then,  the  value  of  vorticity  at  the  wall  needs  be  specified.  However,  there  is 
no  explicit  expression  for  uimniv  —  0)  because,  as  indicated  by  Sutera  (1965), 
the  vorticity  equation  is  coupled,  through  the  no-slip  boundary  condition, 
with  the  velocity  equations.  To  obtain  the  correct  value  of  vorticity  on  the 
wall,  the  equation  for  normal  velocity  Vmn  must  be  solved  first.  By  (6.13) 
and  (6.14),  it  follows 

in  -  ^0  ^mn  =  -iukoUmn  “  with  Umn(O)  =  =  0.  (6.57) 

Using  the  method  of  variation  of  parameters,  the  solution  of  Vmn  satisfying 
the  above  boundary  conditions  is 

f,nhot)  trt  p-nkort  fV 

Vmn  = - ^  J  {i  (^mnAUmn)  dr}'  +  ^  J  (ia;,„„-«„„)  df)' 

“  “  (6.58) 
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But  (6.58)  contains  ximn  which  couples  with  Vmn  through  (6.15)  and  has  not 
been  solved.  Nevertheless,  by  requiring  that  Vmn  is  bounded  as  r;  oo, 
further  analysis  (see  Appendix  A)  shows  that  the  second  boundary  condition 
for  Umn  can  be  expressed  as 


roo 

=  0  (6.59) 

Together  with  (6.56),  we  now  have 


(^mn  +  M2  Dmn  —  1  (6.60a) 

■^1  (^mn  d"  ^2  Dfjin  =  0.  (6.60b) 


Solving  (6.60)  for  C^nn  and  yields 

D  - 

(6.61) 

M2/1-M1/2’ 

where 

Ml  = 

—  .moQ  1 

2A  ■  2  ’  2’  - 

(Wo  - 
2  '' 

(6.62a) 

M2  = 

,,/2A-n^  .man  3 

2A  •  2  '  2’  - 

(6.62b) 

h  = 

.mao  1 

y.  2A  -'—’a 

e-o-d. 

(6.62c) 

h  = 

[°°  »^/2A  — .mao 

y.  2A  -  2  ■ 

3.  ,,  . 

-nfco*7(g^2d) 

Notice  that  in  (6.60)  the  vorticity  at  inflow  has  been  chosen  as  u)mn{Ho)  =  1 
since  the  equation  is  linear.  The  value  of  u)mn  therefore  represents  the  factor 
by  which  the  initial  vorticity  is  amplified  or  attenuated  as  it  approaches  the 
wall.  Before  we  discuss  the  asymptotic  behaviors  of  C^n  and  a  few 
remarks  on  their  significance  and  qualitative  behavior  are  in  order.  First, 
(6.55)  indicates  that  for  a  vortical  disturbance  specified  upstream  by  A  and 
ao,  the  induced  vorticity  on  the  wall  is  directly  related  to  the  value  of  Cmn 
and  Dfnn-  In  the  case  of  A  >•  n^,  jC^nl  represents  the  amplitude  of  vorticity 
at  the  wall,  i.e.  |wTOn(0)|  ~  \Cmn\i  and  represents  the  amplitude  of  nor¬ 
mal  derivative  of  the  vorticity  at  the  wall,  i.e.  K„(0)|  -  \Dmn\.  This  can 
also  be  seen  more  clearly  from  figure  6.10.  Consequently,  the  amplitude  of 
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Cmn  and  J^mn  indicates  the  effectiveness  for  the  disturbance  to  penetrate  the 
boundary  layer  and  modify  the  mean  flow  structures  near  the  wall.  Second, 
the  boundary  condition  (6.59)  shows  that  for  steady  disturbance,  must 
change  sign  in  the  flow  domain;  vorticity  with  sign  opposite  to  the  incom¬ 
ing  vorticity  needs  be  developed  at  the  wall  in  order  to  satisfy  the  no  slip 
boundary  condition.  This  can  be  observed  clearly  form  the  vorticity  contours 
in  figure  6.6.  Note  that  the  amplitude  of  the  induced  vorticity  at  the  wall 
ti'mnCO)}  as  shown  in  figure  6.9,  approaches  a  constant  when  A  becomes  large. 
The  origin  and  implications  of  this  interesting  behavior  of  a;„,„(0)  can  be 
understood  through  the  asymptotic  analysis  for  large  A  and  small  <7o,  and  is 
presented  in  the  next  section. 

6.5  Asymptotic  Behavior 

In  this  section,  the  asymptotic  behavior  of  Cmn  and  Dmn  in  the  vorticity 
expression  is  analyzed  for  the  large  scale  A  1  and  low  frequency  a©  1 
disturbances.  Based  on  this  analysis,  a  scaling  relation  for  the  heat  transfer 
enhancement  in  stagnation  point  flows  in  the  presence  of  upstream  distur¬ 
bances  is  also  derived. 

6.5.1  Vorticity  asymptotes 

The  asymptotic  expression  of  Cmn  for  large  A  and  small  <7o  may  be  obtained 
as  (see  the  Appendix  B) 


2  2 
71^  Of 

\Cmn\  -  lfo[H-(ai-ln^o)y][l--^o^  (6.63) 

where  Oi  =  5(ln2  -1-  7)  and  7  =  0.5772156- ••  is  the  Euler  constant.  The 
expression  for  |jDmn|  can  be  similarly  obtained. 

An  important  observation  for  (6.63)  is  that  the  amplitude  of  the  vortic¬ 
ity  induced  at  the  wall,  up  to  the  leading  order,  is  linearly  dependent  on 
the  normal  distance  Hq  between  the  wall  and  the  upstreani  location  where 
the  disturbance  is  introduced.  This  means  that  for  large  A  and  small  oo, 
within  the  linear  dynamics  regime,  the  amplification  factor  due  to  the  vortex 
stretching  for  the  initial  vortical  disturbance  has  an  upper  limit  set  by  Hq. 
This  explains  the  behavior  seen  in  figure  6.9  as  A  becomes  large.  FVom  (6.63), 
one  also  finds  that  the  unsteadiness  of  the  disturbances  tends  to  decrease  the 
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induced  vorticity  amplitude  at  the  wall  compared  with  the  steady  case,  but 
only  to  the  second  order  in  terms  of  frequency  oo-  In  figure  6.13,  the  vorticity 
profiles  at  different  frequencies  are  computed  numerically  for  A  =  4  to  show 
the  effect  of  the  disturbance  unsteadiness.  Notice  that  although  the  numeri¬ 
cal  computation  uses  the  full  nonlinear  equations,  the  vorticity  value  at  the 
wall  indeed  follows  the  asymptotic  behavior  (6.63)  and  shows  only  modest 
changes  relative  to  steady  case  for  oq  <  0.5.  FVom  a  viewpoint  of  time  scale, 
oq  <  1  implies  that  the  disturbance  turn-over  time  is  much  longer  than  the 
time  scale  of  mean  flow  distortion;  vortex  stretching  is  predominant,  and 
hence  generating  flow  structures  similar  to  the  steady  case.  For  high  fre¬ 
quency  case  00  1,  (6.28)  is  not  valid,  and  in  the  course  of  disturbances 

being  convected  towards  the  wall,  many  cycles  of  oscillation  have  completed. 
The  net  vorticity  induced  at  the  wall  is  small  due  to  the  cancellation  effect 
of  the  incoming  disturbance  vorticity  with  alternating  signs.  This  effect,  also 
known  as  ‘vortex  piling’,  has  been  analyzed  by  Hunt  (1973).  Our  numeri¬ 
cal  results  in  figure  6.13  shows  that  indeed  the  induced  vorticity  at  the  wall 
decays  monotonically  with  <to  and  becomes  rather  small  when  oq  >  2.5. 

6.5.2  Heat  transfer  scaling 

With  the  amplitude  of  the  velocity  disturbance  at  the  inflow  Vm„{Ho)  kept 
constant,  the  vorticity  Umn  introduced  at  the  inflow  varies  as  a  function  of  A 

^mnifJo,X)  —  U}mn{Ho,l)  X  *  (6.64) 

where  Umn{Ho,  1)  is  a  reference  value  of  vorticity  for  A  =  1.  Recall  that,  Cmn, 
Dmm  functions  of  A  and  oo,  represent  the  amplification  factor  for  the  wall 
vorticity  value.  The  amplification  factor  increases  with  A,  initially  at  a  large 
rate  but  eventually  approaches  a  constant  value  determined  by  Hq.  Together 

with  the  A  *  factor,  this  gives  rise  to  an  overall  optimum  amplification  for 
A  4 

Once  the  vorticity  is  obtained,  the  corresponding  velocity  disturbance  can 
be  found  by  (6.58)  as 

^-nkoT)  pt 

Vmn  =  {iUmn  -  u^n)  dr}.  (6.65) 

As  shown  in  appendix  A,  the  second  term  of  the  integrand  in  (6.65)  is  ex¬ 
pected  to  be  much  smaller  than  the  first  for  Aiq  <C  1.  Hence  the  dependence 
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of  Vmn  on  A  directly  follows  that  of  Umn-  When  A  is  large,  the  vorticity  am¬ 
plification  factor  Cmn  approaches  a  constant.  Accordingly,  |nlmox  A“^/^,  a 
trend  noted  in  figure  6.8.  Now  from  (6.23),  it  follows  that 


Ah  _  ^o(O) 
h  ~  0'{O) 


^  ^  l^mnl  l^mnl 

m,n 


(6.66) 


For  the  temperature  disturbance  Omni  the  linearized  equation  of  (6.17)  is 

+  Pr<j)  -  {n^kl  +  im  Prao)Omn  =  PrVmi0  (6-67) 


The  amplitude  of  the  temperature  fluctuation  6mn  can  be  deduced;  viz.  dmn 
will  be  proportional  to  the  amplitude  Vmn  and  follows  the  same  asymptotic 
dependence  on  A  as  u>mn  and  v^n-  Hence,  we  have 

~  ~  ~  (6.68) 

For  the  low  frequency  case  oo  <IC  1,  the  heat  transfer  enhancement  scales 
with  the  disturbance  parameters, 

X  ~  (6-69) 

m,n 

When  A  :»  1,  we  have  Ah/h  ~  A“^  for  fixed  Ap,  a  behavior  noted  in  figure 
6.8.  We  can  rewrite  ^  in  terms  of  the  disturbance  length  scale  1^,  and  obtain 
the  following  scaling  relation  for  the  relative  heat  transfer  enhancement  at 
low  disturbance  frequency 


Ah 

h 


5 


|Vmn(^fo)P 


m,n 


(6.70) 


The  above  expression  reveals  different  roles  played  by  the  various  flow  param¬ 
eters  in  heat  transfer  enhancement.  First,  the  enhancement  is  proportional 
to  the  square  of  the  disturbance  amplitude  Ap  due  to  the  net  convective  flux 
by  the  disturbance  modes.  Second,  the  length  scale  of  the  disturbance  has 
a  critical  effect  upon  the  heat  transfer  enhancement.  For  large  scale  distur¬ 
bances,  the  heat  transfer  enhancement  decreases  with  increased  length  scale. 
The  most  effective  disturbance  will  be  those  with  length  scales  comparable 
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to  the  boundary  layer  thickness.  Finally,  the  overall  effect  of  the  frequency 
of  the  disturbances  is  to  reduce  the  heat  transfer  enhancement.  The  decrease 
is  of  the  second  order  when  the  frequency  is  low,  and  at  high  frequencies  the 
heat  transfer  enhancement  approaches  to  zero  owing  to  the  rapid  decay  of 
the  incoming  disturbances.  So  far  the  Prandtl  number  has  been  kept  at  a 
constant  in  the  analysis  Pr  =  0.71.  For  a  different  Prandtl  number,  both  the 
mean  thermal  boundary  layer  thickness  and  the  amplitude  of  the  will 
be  affected.  Nevertheless,  the  qualitative  dependence  on  Prandtl  number  can 
be  inferred  from  (6.23).  For  small  Pr,  the  Ahfh  increases  with  the  increase  of 
Pr,  but  as  Pr  becomes  very  large,  the  effect  diminishes  due  to  the  extremely 
thin  thermal  boundary  layer.  So  there  will  be  a  optimal  value  of  Prandtl 
number  for  each  fixed  value  of  X  which  gives  the  maximum  heat  transfer  en¬ 
hancement.  Numerical  calculations  show  that  this  optimal  Prandtl  number 
decreases  with  increasing  A  from  Pr=  7  at  A  =  1.5  to  Pr  =  0.71  at  A  =  12. 
However,  for  a  fixed  Prandtl  number,  the  optimal  value  of  A  is  still  around 
A  =  4,  and  the  overall  maximum  value  of  Ah/h  occurs  around  A  =  4  and 
Pr=  1.5  (Xiong,  2004). 

6.6  Discussion  of  Pree-Stream  Turbulence 

One  of  the  primary  goals  of  the  preceding  analysis  is  to  provide  an  improved 
understanding  on  the  free-stream  turbulence  effect  in  stagnation  point  flows. 
In  this  section,  we  extend  the  incoming  disturbances  from  discrete  harmonic 
modes  to  free-stream  turbulence.  The  latter  is  assumed  to  be  isotropic  and 
homogeneous.  This  extension  allows  a  comparison  to  experimental  measure¬ 
ments.  This  also  serves  as  a  test  of  applicability  of  the  analysis  to  relevant 
engineering  problems.  The  formulation  is  analogous  to  that  in  RDT,  i.e. 
the  overall  changes  of  the  turbulence  statistics  are  obtained  by  integrating 
over  all  the  Fourier  modes  once  the  modal  distortion  for  each  of  them  is 
known.  The  goal  is  to  derive  a  scaling  relation  between  the  characteristics 
of  the  free-stream  turbulence  and  the  relative  heat  transfer  enhancement  at 
the  wall. 

On  expressing  the  velocity  fluctuations  as 

//£  U  (/Cl,  /C2,  /C3)  exp  [i(/Ci^  -f  K2T}  -|-  -I-  at)]dKidK2dK3 

(6.71) 
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and  using  Taylor’s  hypothesis,  the  free-stream  turbulence  is  treated  as  ‘frozen’ 
while  being  convected  into  the  domain  through  the  inflow  boundary  with  a 
time  frequency  a  =  -Uk2-  Let  r}  =  Hohe  the  location  of  the  inflow  bound¬ 
ary,  then 


n((,  Ho-m,0  =  j  f  r  K2,  k,)  dn,. 

(6.72) 

While  being  convected  towards  the  stagnation  point,  the  free-stream  tur¬ 
bulence  experiences  an  accumulated  vortex  stretching  in  ^  direction  by  the 
diverging  mean  flow.  As  a  result,  this  leads  to  ki  <C  k2,  k^  in  the  stagnation 
point  region.  Thus  the  dependence  of  u  on  ^  may  be  neglected  in  the  above 
expression,  and  the  inverse  transformation  of  Ui„  is  now  written  as 

u*,(-C2,  Ks)  =  J jT  -  Ut,  C)  d(dt  (6.73) 

The  heat  transfer  enhancement  in  the  presence  of  free-stream  turbulence  can 
be  computed  in  terms  of  the  downstream  velocity  and  temperature  spectra 
n{r),K2,K3),  0{r],K2,K3).  Following  the  discussion  in  preceding  sections,  we 
assume  that  for  the  77-velocity  v  and  temperature  0 

v{t),  K2,  K3)  =  G{ri)v{Ho,  K2,  K3)  ~  0{t},  K2,  ks)  (6-74) 

where  function  G{t})  represents  the  amplification  ratio  at  downstream  loca¬ 
tion  77  for  mode  v{k2,  K3)  due  to  mean  flow  straining  and  viscous  dissipation. 
Taking  the  amplitude  of  G{ti)  as  C{k2,K3),  i.e. 

\G{r),K2,K3)\  ~  |C'(k2>^3)|  (6.75) 

The  overall  effect  of  free-stream  turbulence  on  the  heat  transfer  can  be  es¬ 
timated  by  summing  the  contribution  from  all  the  modes  of  different  wave 
numbers  and  frequencies.  In  view  of  (6.23),  for  a  single  mode  disturbance 
specified  at  the  inflow,  the  heat  transfer  enhancement  mainly  results  from 
the  second  order  interaction  terms.  Thus  in  the  present  linear  analysis  when 
(6.23)  and  (6.74)  are  used  to  find  the  contribution  from  mode  Vmn,  the  non¬ 
linear  term  A/4  can  be  written  as 

M^i{K2,tZ3)  ~  lC’(^2j  ^3)P  |u(«2>  ^3)1^  |C'(«2>  ^3)1^  ^n(^2>  ^3)  (6.76) 
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where  the  $«(/C2,  «3)  is  the  energy  spectra  density  for  the  free  stream  tur¬ 
bulence.  Since  the  free-stream  turbulence  is  assumed  to  be  isotropic  and 
homogeneous,  the  energy  spectrum  E{k)  is  defined  as: 

E{k)  =  27r$«(«)/c2  77^ 

The  turbulence  intensity  and  the  integral  length  scale  L  may  be  related 
to  the  amplitude  and  fundamental  wavenumber  of  energy  containing  range 
by 

Al  ~  ~  ^  (6.78) 

The  total  contribution  from  all  wave  number  components  is  the  integration 
of  the  turbulent  energy  spectrum,  i.e. 

^pI]|«(^o,«2,/C3)P  ~  ///  ^ii{K)(^K  -  Tj  r^dinL) 

K2yK3  J  J  Jo  Jo  q^L  ' 

(6  T9) 

Substitute  (6.76),  (6.79)  into  the  expression  for  heat  transfer  enhancement 
(6.23),  it  becomes 


~  I?  f  E{^L)\C(kL)\U{KL)  (6.80) 

where  E{kL)  =  E{k)/ (q^L).  In  most  of  the  engineering  problems  involving 
turbulence  impinging  on  a  stagnation  point,  the  turbulent  eddy  turn-over 
time  is  t3T)ically  much  longer  than  the  time  scale  associated  with  the  mean 
flow  straining,  i.e.  a©  <  1.  In  addition,  the  turbulence  length  scale  is 
assumed  much  larger  than  the  boundary  layer  thickness,  i.e.  L  >  5.  By 
(6.68),  for  a  fixed  upstream  location  it  follows 

\C{kL)\  ^  ^  {kL){£)  for  0  <  K  <  Kmax  (6.81) 

Hence, 

Ahr 

—J—  ^  /  E{kL)  {nLf  d{KL)  (6.82) 

Jo 

where  k^ax  represents  the  highest  wavenumber  component  having  contribu¬ 
tion  to  the  heat  transfer  augmentation.  From  the  previous  analysis,  the  heat 
transfer  enhancement  reaches  its  maximum  value  with  the  disturbance  scale 
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similar  to  the  the  boundary  layer  thickness.  At  even  smaller  scales  the  dis¬ 
turbances  decay  rapidly  because  of  the  viscous  dissipation  and  have  little 
effect  on  the  wall  quantities. 

For  practical  engineering  problems,  such  as  flows  over  a  bluff  body,  the 
boundary  layer  thickness  at  the  stagnation  region  scales 


Iq  =  y/l/*/A*  ~  ~ 


D 

■/fle 


(6.83) 


where  the  Reynolds  number  is  based  on  the  mean  flow  U  and  the  diameter  of 
curvature  D,  see  figure  6.1a.  For  an  effective  strain  rate  A*,  e.g.  A*  =  4U/D 
for  a  circular  cylinder,  Kmax  roay  be  expressed  as 

KmaxL  ~  ^  ~  (6.84) 

Furthermore,  the  energy  spectrum  of  the  free  stream  turbulence  is  as¬ 
sumed  to  follow  the  Kolmogoroff  —5/3  law 


E{kL) 


(6.85) 


Substituting  (6.84)  and  (6.85)  into  (6.82),  the  heat  transfer  enhancement  at 
stagnation  point  can  now  be  correlated  to  the  free-stream  turbulence  param¬ 
eters  by 


h 


(6.86) 


Note  that  here  the  three-dimensional  spectrum  function  is  used.  Following 
the  form  of  organized  disturbances  in  previous  discussion,  a  one-dimensional 
spectra  is  probably  more  appropriate.  However,  at  high  Reynolds  number, 
it  also  follows  the  same  —5/3  law  as  the  three  dimensional  spectrum. 

In  order  to  examine  heat  transfer  scaling,  (6.86)  is  compared  against 
the  experimental  measurements  for  stagnation  point  flows  in  the  presence  of 
free-stream  turbulence.  A  recent  experiment  was  conducted  by  Ames,  Wang, 
and  Barbot  (2002),  in  which  the  heat  transfer  to  a  model  vane  is  measured 
for  six  different  inlet  turbulence  conditions  with  turbulence  intensity  up  to 
14  percent.  The  different  characteristics  of  the  free-stream  turbulence  are 
generated  using  mesh  biplanar  grid  and  various  mock  combustion  system 
configurations.  The  experimental  set  up  is  representative  to  modern  dry  low 
NOx  and  aeroderivative  combustors.  The  downstream  vane  heat  transfer 
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measurements  can  serve  as  a  database  for  the  validation  of  predictive  meth¬ 
ods.  Another  reason  to  choose  this  experiment  is  that  the  characteristics  of 
free-stream  turbulence  i.e.  intensity  and  length  scale  are  measured  with  the 
model  vane  present  in  the  flow.  Thus  the  mean  velocity  follows  the  Hiemenz 
profile  at  the  locations  where  the  turbulence  is  measured.  This  is  the  same 
case  as  has  been  assumed  in  the  present  analysis.  In  some  other  experiments, 
e.g.  Van  Fossen  et  al.  (1995),  the  free  stream  turbulence  is  first  measured 
m  the  uniform  flow  without  the  model.  When  the  model  is  present,  the  tur¬ 
bulence  characteristics  are  obtained  by  extrapolation  using  a  power  law  for 
decaying  turbulence.  At  locations  close  to  the  stagnation  point,  the  turbu¬ 
lence  characteristics  obtained  by  this  method  will  be  significantly  different 
from  those  obtained  with  the  presence  of  the  model.  Although  the  differ¬ 
ence  becomes  smaller  beyond  a  distance  on  the  order  of  D  away  from  the 
wall^  the  Hiemenz  flow,  as  a  good  approximation  at  the  stagnation  point  for 
bluff-body  flows,  is  also  only  valid  within  the  order  of  D  away  from  the  wall. 
Hence,  in  figure  6.14  the  present  heat  transfer  correlation  is  compared  with 
the  measurement  of  Ames  et  al.  (2002)  at  four  high  turbulence  levels  gener¬ 
ated  by  grid  as  well  as  by  aero-derivative  and  dry  low  NOx  mock  combustor 
system.  The  Reynolds  number  based  on  the  leading  edge  diameter  ranges 
from  58,000  -  232,000;  turbulence  integral  scale  from  O.llD-  l.OD  and  the 
turbulence  intensity  from  8  -  14%.  The  collective  experimental  uncertainty 
IS  ±5  percent  for  turbulence  measurement  and  ±3  percent  for  heat  trans¬ 
fer  data.  Although  some  scatter  is  present,  the  correlation  appears  to  agree 
reasonably  well  with  the  experimental  data  for  the  correlation  parameter  H 
over  the  range  5  to  35.  Notice  that  the  turbulence  levels  in  these  experiments 
are  quite  high,  but  the  analysis  indicates  that  only  small  scale  components 
contribute  effectively  to  the  heat  transfer  enhancement.  So  even  the  total 
turbulence  level  is  high,  the  amplitude  at  the  small  scale,  i.e.  the  scale  of 
boundary  layer  thickness  that  affects  the  heat  transfer  most,  would  still  be 
relatively  small.  For  a  single  mode  velocity  disturbance  at  a  level  equivalent 
to  10  percent  of  the  mean  flow,  the  numerical  results  show  that  the  dis¬ 
turbance  evolution  can  still  be  largely  described  by  linear  vortex  dynamics. 

his  may  explain  why  the  correlation  based  on  linear  analysis  seems  to  hold 
even  for  the  case  of  high  turbulence  intensity.  It  is  also  interesting  to  notice 
that  the  present  correlating  parameter  is  close  to  the  square  of  the  empir¬ 
ical  TLR  parameter  TLR  =  T^Re%^\D/L^y/^  proposed  by  Ames  (1997), 
if  the  turbulence  integral  length  scale  were  replaced  by  the  ‘energy  scale’ 
Lu  =  1.5|ti'l7cT,  where  |u'|  is  the  r.m.s.  streamwise  fluctuation  velocity  and 
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cr  is  the  turbulent  dissipation  rate. 
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Figure  6.1;  a)  Disturbed  stagnation  point  flow  at  the  leading  edge  of  a  two- 
dimensional  bluff  body,  b)  Hiemenz  boundary  layer  flow  with  upstream  in¬ 
coming  disturbances. 


Figure  6.3:  Streamline  patterns  in  th 
A.  Flow  is  downward  and  the  spanwi 
A  =  1.1, 2, 4, 16,32,64  from  a  to  f. 


stagnation  plane  for  different  values  of 
i  width  is  one  disturbance  wavelength. 
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Figure  6.4:  Streamline  patterns  in  the  stagnation  plane  for  different  values 
of  i4p  at  A  =  4.  Flow  is  downward  and  the  spanwise  width  is  one  disturbance 
wavelength.  Ap  =  0.02, 0.05, 0.07, 0.10, 0.12, 0.15  from  a  to  f. 


176 


T1 


Figure  6.5:  Profiles  for  steady  {a  =  0)  velocity  and  temperature  disturbances 
at  A  -  4.  Solid  line  is  the  fundamental  mode,  The  dash,  dash  dot  and  dotted 
line  correspond  to  mean  flow  modification,  second  and  third  harmonics  for 
«mn,  Vmn,  Bmn,  and  second,  third  and  fourth  harmonics  for  Umn  and  Wmn- 
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Figure  6.6:  Vorticity  contours  in  stagnation  plane  at  A  =  4,  Flow  is  downward 
and  the  spanwise  width  is  27r.  Solid  lines  are  for  the  positive  values  of  co  and 
dash  lines  for  the  negative. 
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Figure  6.7:  Temperature  contours  in  stagnation  plane  at  A  =  4.  Flow  is 
downward  and  the  spanwise  width  is  27r. 
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Figure  6.9:  Vorticity  magnitude  at  the  wall 
fundamental  mode. 
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Figure  6.10:  A  typical  vorticity  profile  composed  of  two  confluent  hyper¬ 
geometric  functions.  A  =  4,  oo  =  0 
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Figure  6.11:  Comparison  of  vorticity  profile  between  numerical  and  asymp¬ 
totic  solutions  at  A  =  4. 
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Figure  6.12:  Comparison  of  vorticity  profile  between  numerical  and  asymp¬ 
totic  solutions.  A  =  36. 
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Figure  6.14:  Experimental  data  of  leading  edge  heat  transfer  enhancement 
under  free-stream  turbulence  correlated  by  the  present  scaling  parameter  %. 


Chapter  7 

Summary  and  Conclusion 


We  have  developed  a  numerical  method  for  three-dimensional  compressible 
flow  based  on  fourth  order  finite  difference  scheme  in  curvilinear  coordinates, 
and  a  very  efficient  dual-time  linearized  subiteration  scheme  for  fully  implicit 
time  marching.  This  method  is  suitable  for  numerically  simulating  turbulent 
flows  around  bluff  bodies.  It  has  been  implemented,  validated,  and  used  to 
perform  large  eddy  simulations  of  free-stream  turbulence  impinging  on  an 
elliptical  leading  edge  and  the  resulting  heat  transfer  enhancement. 

A  new  blending  procedure  is  developed  by  which  independent  but  sta¬ 
tistically  identical  realizations  of  homogeneous  isotropic  turbulence  are  com¬ 
bined  into  a  unified  field  to  represent  free-stream  turbulence  realistically. 
The  blending  is  performed  over  a  blending  zone  at  the  interface  between  dif¬ 
ferent  realizations  using  a  smooth  varying  function.  Extra  dilatation  field 
within  the  blending  zone  is  further  removed  using  Helmholtz  decomposition 
theorem.  Outside  the  blending  zone,  the  original  turbulence  fields  remain 
unchanged.  By  construction,  the  combined  field  preserves  the  turbulence 
intensity,  and  the  change  to  other  statistical  quantities  are  shown  to  be  min¬ 
imal.  The  method  is  simple  yet  effective  and  could  be  useful  in  other  direct 
or  large  eddy  simulations  in  which  effects  of  sustained  free-steam  turbulence 
are  important. 

The  large  eddy  simulation  results  characterize  the  anisotropy  of  turbu¬ 
lence  due  to  the  mean  flow  strain  as  it  approaches  the  stagnation  point.  The 
turbulent  budgets  for  the  Reynolds  stresses  in  such  a  strongly  strained  flow 
are  presented,  and  different  roles  in  the  dynamics  for  different  terms  are  dis¬ 
cussed.  These  results  provide  guidelines  for  modelling  strongly  strained  tur¬ 
bulent  flows  for  Reynolds  averaged  Navier-Stokes  equations,  a  subject  which 
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has  been  proved  difficult  in  the  past.  The  numerical  results  also  shows  that 
free-stream  turbulence  significantly  elevates  the  surface  heat  transfer  coeffi¬ 
cients,  but  the  distribution  of  the  heat  transfer  coefficient  along  the  surface 
IS  essentially  kept  the  same  as  in  the  laminar  case.  Both  of  these  results 
Jow  good  agreement  with  the  corresponding  experimental  measurements. 

he  numerical  simulation  also  reveals  that  small  scale,  intense  vortical  flow 
structures  generated  at  the  leading  edge  by  vortex  stretching  induces  signif¬ 
icant  changes  in  the  thermal  boundary  layer  thickness,  and  are  the  primary 
cause  of  the  observed  heat  transfer  enhancement. 

In  the  theoretical  study,  we  analyze  the  distortion  of  unsteady  three  di¬ 
mensional  organized  disturbances  in  an  incompressible  Hiemenz  boundary 
layer  and  its  effect  on  the  wall  heat  transfer.  Based  on  linearized  disturbance 
equations,  it  is  found  that  the  vorticity  outside  the  boundary  layer  can  be 
expressed  analytically  in  terms  of  confluent  hyper-geometrical  functions,  pa¬ 
rameterized  by  the  disturbance  length  scale  and  temporal  frequency.  When 
the  scale  of  the  disturbance  is  large  and  the  frequency  is  low,  an  approximate 
^ymptotic  solution  is  obtained  with  explicit  dependence  on  the  disturbance 
length  scale  and  frequency.  This  solution  compares  well  with  the  full  nonlin¬ 
ear  numerical  solutions  over  a  wide  range  of  disturbance  parameters. 

It  IS  further  shown  that  the  ratio  between  the  disturbance  length  scale 
and  the  boundary  layer  thickness  is  the  critical  parameter  in  determining 
the  amplification  factor  of  incoming  vorticity.  Essentially,  it  represents  the 
interaction  between  the  vortex  stretching  effect  and  viscous  diffusion.  The 
amplification  factor  is  found  to  be  inversely  proportional  to  the  length  scale 
except  at  very  small  scales  where  it  increases  with  increased  length  scale. 
The  maximum  amplification  is  found  for  the  disturbances  with  a  length  scale 
about  five  times  the  Hiemenz  boundary  layer  thickness.  The  associated  heat 
transfer  enhancement  also  strongly  depends  on  the  disturbance  length  scale 
and  IS  analyzed  through  the  induced  vorticity  at  the  wall.  Compared  to  the 
steady  ca^,  the  heat  transfer  enhancement  is  reduced  by  the  unsteadiness  of 
the  disturbance  but  the  effect  is  of  second  order  when  the  frequency  is  low. 

The  analysis  is  further  extended  to  the  case  of  homogeneous  and  isotropic 
free-stream  turbulence.  The  turbulence  energy  spectrum  is  assumed  to  follow 
the  Kolmogoroff’s  -5/3  law  and  the  integral  scale  is  much  larger  than  the 
boundary  layer  thickness.  Under  these  conditions,  a  new  scaling  correlation 
IS  derived  between  the  relative  heat  transfer  enhancement  and  the  turbu¬ 
lence  intensity,  integral  length  scale  and  the  mean  flow  Reynolds  number. 

In  comparison  to  recent  experimental  data  on  turbine  blade  heat  transfer 
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in  the  presence  of  free-stream  turbulence,  the  present  correlation  provides  a 
reasonable  guide  to  the  observed  variations. 

This  technical  report  represents  a  status  report  at  the  time  of  writing. 
The  post  processing  of  the  simulation  data  is  continuing  and  the  finalized 
version  including  all  the  results  and  discussions  will  be  reported  in  the  PhD 
thesis  of  Xiong  (2004). 
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Appendix  A 

Incompressible  Potential  Flow 


We  derive  the  incompressible  two  dimensional  potential  flow  solution  of  a 
point  source  placed  midway  between  two  parallel  planes  with  a  uniform  in¬ 
coming  stream,  see  figure  C.l.  The  configuration  is  prototypical  to  experi¬ 
ments  where  a  blunt  leading  edge  model  with  flat  plate  after-body  placed  in 
a  wind  tunnel  test  section,  and  due  to  the  thickness  of  the  model  in  com¬ 
parison  with  the  height  of  the  test  section,  the  flow  passage  is  significantly 
blocked.  The  blockage  causes  the  flow  field  to  differ,  particularly  in  pressure 
distribution,  from  the  open  flow  case  where  the  walls  are  absent.  To  make  a 
direct  comparison  between  computation  and  experiment,  the  blockage  effect 
often  needs  to  be  accounted  explicitly  in  the  computational  configuration. 

Let  Uoo  be  the  uniform  incoming  velocity  at  far  upstream,  a  be  the  half 
height  of  the  test  section,  and  6  the  half  thickness  of  the  model.  The  coordi¬ 
nate  is  chosen  such  that  the  forward  stagnation  point  is  at  the  origin  x  =  0. 
The  goal  here  is  to  find  the  strength  27rm  and  the  location  x  =  c  of  a  two 
dimensional  point  source  S  so  that  at  far  downstream,  the  required  blockage 
ratio  hfa  can  be  reproduced. 

Following  Milne-Thomson  (1960),  the  complex  potential  for  a  point  source 
placed  midway  between  two  planes  with  a  general  uniform  incoming  stream 
U  can  be  written  as 

W  =  (f)  +  iip  =  —Uz  —  m  In  sinh  — —  (A.l) 

2a 

here  (f>  and  ^  are  the  velocity  potential  and  stream  function,  and  z  =  x  +  iy. 
By  the  definitions  of  complex  function  In  z  and  sinh  z 

\nz  =  ln(\/x^  +  y^)  +i  arctan  —  (A.2) 

X 
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sinhz  =  cosysinha:  +  i  sin  ycoshx 
The  stream  function  in  (A.l)  is 

Ip  =  -Uy  -  m arctan  [tan  ^  coth  — ^1 

2a  2a  ' 

The  dividing  streamline,  which  forms  the  boundary  of  the  Rankine  body, 
corresponds  to  =  0.  This  requires 


(A.3) 


(A.4) 


.  Uy  ^  Try  ,  7r(x  —  c) 

tan  —  =  —  tan  —  coth  — ^ - - 

m  2a  2a 


At  far  downstream  x  -xx),  coth  1,  thus 


whose  solution  is 


^  Uy  Try 

tan  —  =  —  tan  — 
m  2a 


Uy  Try  . 

—  =  -■ ^inTT,  n  =  0,1,2... 
m  2a 


(A.5) 


(A.6) 


(A.7) 


Physical  conditions  indicate  that  only  n  =  0  and  n  =  1  are  needed,  so  that 
at  far  downstream  the  y  coordinates  of  the  Rankine  body  are 


y  =  0;  y  =  ±6  =  db 


2am7r 


as  X  — >  oo 


(A.8) 


2aU  4-  myr 

Note  that  to  ensure  |y|  <  a  at  x  — >  oo,  the  incoming  velocity  U  must  satisfy 

Trm 


U  > 


2a 


Furthermore,  the  streamwise  velocity  u  can  be  obtained  from  xp  by 


coth^iif::^! 

2tt 


—  u 

dy  ~  ~^cos^^)  +  sin2(M)  coth^ 


Setting  «(0, 0)  —  0,  we  obtain  the  location  of  the  source  c 

a  ,  ,  2aU  +  Trm , 


But,  if  we  now  let  x 


- ) 

TT  '2af7  — Trm'' 

-oo, (A.IO)  leads 

TT  y  rr 

u^V-—^V 


(A.9) 


(A.10) 


(A.11) 


(A.12) 
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That  is,  the  resulting  velocity  at  upstream  differs  from  the  value  U  which 
we  initially  imposed!  This  apparent  contradiction  stems  from  the  fact  that 
due  to  the  presence  of  the  planes,  the  velocity  generated  from  the  source, 
which  would  vanish  asl/ratr->ooin  open  flows,  does  not  vanish  in  the 
present  case;  instead,  it  reaches  a  constant  value  at  x  =  ±oo.  Indeed,  by 
continuity  equation  and  in  the  absence  of  imposed  mean  flow  U,  a  point 
source  of  strength  27rm  placed  between  two  planes  y  =  rto  generates  velocity 

u  =  ±-—  at  X  =  ±oo  (A.13) 

2a 

at  X  =  ±oo.  The  incoming  velocity  at  x  =  — oo  is  therefore  a  superposition 
of  the  imposed  incoming  mean  flow  U  and  the  outgoing  velocity  ^  from 
the  source  inside  the  domain.  If  the  boundary  condition  is  set  u  =  Uoo  at 
X  =  —  oo,  then  we  need  choose 

U  =  U^  +  ^  (A.14) 

for  the  resulting  velocity  to  be  compatible.  Hence,  for  given  lengthes  a, 
b,  and  velocity  Uoo,  substituting  (A.14)  into  (A.8)  and  (A. 11),  we  find  the 
required  downstream  blockage  can  be  realized  by  placing  the  point  source  S 


of  strength  27rm,  where 

ab 

/  u\^oo 

7r(a  —  b) 

(A.  15) 

at  location  x  =  c,  where 

a  ,  'Km  . 

TT  Oi^oo 

(A.16) 

An  example  of  the  streamline  pattern  for  a  =  5,  6  =  1  and  Uoo  =  1  is  shown 
in  figure  A.  2.  Note  that  obviously  in  the  potential  framework,  the  boundary 
layers  developed  along  the  model  surface  and  the  wind  tunnel  walls  are  ne- 

glected.  Furthermore,  in  the  near  field,  the  generated  Rankine  body  shape 
does  not,  in  general,  follow  the  specific  geometry  of  the  leading  edge  of  a  test 
model.  The  solution,  however,  can  serve  as  an  initial  condition  for  numerical 
computations  where  specific  leading  edge  geometry  is  incorporated,  and,  in 
low  Mach  number  case,  provides  a  reasonably  good  boundary  condition  at 
far  field  for  the  compressible  flow  in  question. 
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Figure  A.l:  A  source  with  uniform  incoming  flow  placed  between  two  parallel 
plates 


Figure  A.2:  Streamline  pattern  of  a  point  source  with  uniform  stream  placed 
between  two  planes,  a  =  5,  6  =  1  and  f/go  =  1 


Appendix  B 

On  the  Coefficient  of  Numerical 
Dissipation 


Center  differencing  scheme  (differencing  scheme  on  a  symmetric  stencil)  is 
dispersive  but  not  dissipative.  Owing  to  the  stencil  symmetry,  the  dissipation 
terms  in  the  truncation  error  —  the  even  order  derivatives  are  completely 
canceled  out.  The  discretization  error  in  this  case  consists  of  only  dispersion 
terms  associated  with  the  odd  derivative  terms.  For  a  simple  traveling  wave 
computed  using  center  differencing  scheme,  for  instance,  the  wave  will  travel 
at  a  different  speed  due  to  the  numerical  dispersion,  but  the  amplitude  of 
the  wave  is  maintained. 

The  inherent  non-dissipative  nature  makes  the  center  differencing  scheme 
susceptible  to  numerical  instability,  particularly  in  the  simulations  of  turbu¬ 
lence.  Once  some  high  wavenumber  components  are  generated  and  can  not 
be  represented  adequately  by  the  underlying  grid,  they  tend  to  persist  in  the 
course  of  computation  and  often  eventually  leads  to  numerical  divergence.  A 
common  cure  to  this  problem  is  to  add  controlled  artificial  dissipation  terms 
to  the  governing  equations  to  suppress  the  unwanted  high  wavenumber  com¬ 
ponents.  Another  way  of  dealing  it  is  to  use  numerical  filtering,  which  will 
not  be  discussed  here. 

Specifically,  consider  a  general  one  dimensional  equation  in  its  conserva¬ 
tive  form 

ft  +  FM  =  0  (B.l) 

where  the  subscripts  denote  temporal  and  spatial  partial  derivatives.  Suppose 
in  computational  domain  [0, 1],  a  uniform  grid  consists  of  N  points  is  used. 
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then  Xi  _  ih,  \<i^N,  where  h  =  l/N  is  the  grid  spacing.  A  five  point 
fourth  order  center  differencing  scheme 


~  /»42  —  8(/i_i  — 


+  0{h^) 


(B.2) 


is  chosen  to  approximate  the  spatial  derivative.  To  enhance  the  robustness 

o  di^fference  equation  resulting  from  (B.2),  a  numerical  dissipation  term  of 
the  form 

-‘''■*(0)  (B.3) 

may  be  to  the  right  hand  side  of  (B.l).  Here  o  >  0  is  the  coefficient 
controlling  the  amount  of  numerical  dissipation  added,  whereas  the  h*  keeps 
t  e  dissipation  term,  if  o  is  not  excessively  large,  of  the  same  order  of  the 
trancation  error  in  {B.2).  In  this  way,  it  is  expected  that  the  formal  accuracy 
of  the  difference  equation  will  not  be  ruined  by  the  addition  of  (B.3).  The 

fourth  order  derivative  in  (B.3)  may  be  approximated  by  the  following  seven 
point  4th  order  scheme 

^  ^  ^f“^-^'-3  +  /<+3)  +  12(/,_2i+/j+2)-39(/,_i  +  /<+i)  +  56/J. 

Now  the  problem  reduces  to  how  to  the  specify  the  value  of  a  so  thafthe 
unresolved  scales  can  be  effectively  damped  while  the  resolved  scales  remain 
essentially  unaffected.  In  practice,  the  value  of  a  is  often  chosen,  to  a  large 
extent,  empirically  depending  on  the  computations  in  question.  Here,  we 
analyze  the  proper  bound  of  a  using  modified  wavenumber  analysis* '  the 

purpose  is  to  provide  a  guideline  in  determining  the  value  of  a  in  practical 
computations. 

funlbn^^'"  modified  wave  number  analysis,  let  /  be  a  simple  periodic 

/  =  c’*®  then  fx  =  (B.5) 

The  relative  error  contained  in  the  numerical  differentiation  by  (B.2)  is 


(B.6) 


\fx-fx\  .  16sin(*:/i)  -  2sin(2A:/i) 

HIT  =  * - r2kh  (B-6) 

On  the  other  hand,  the  relative  error  resulting  from  (B.3)  can  also  be  similarly 
expressed  as 


16sin{kh)  -  2s\n{2kh) 
12kh 


-  ^^-2cos(3fch)  +  24cos{2kh)  -  78cos(fch)  +  56 
l/xl  6kh  ^ 


(B.7) 


196 


Figure  B.l:  Relative  error  of  4th  order  center  differencing  scheme  and  4th 
order  dissipation  term  (scaled)  as  a  function  of  hk. 

To  maintain  the  formal  accuracy  of  the  numerical  solution,  we  require  that 
the  maximum  of  the  (B.7)  is  one  order  of  magnitude  smaller  than  maximum 
of  (B.6),i.e. 

max  ah'* I I  <  0.1  max\fx  —  fx\  (B.8) 

For  the  particular  choice  of  difference  schemes  used  here,  it  is  found  that  this 
can  be  achieved  by  setting 

ah  <  0.01  (B.9) 

For  other  differencing  schemes,  similar  results  can  also  be  found.  The  com¬ 
parison  of  (B.6)  and  (B.7)  resulting  from  (B.9)  is  shown  in  the  figure  B.l  as  a 
function  of  kh.  As  can  be  seen,  the  error  in  this  case  caused  by  the  numerical 
dissipation  term  is  negligible  compared  to  the  inherent  difference  truncation 
error.  If  we  take  (B.9)  as  an  upper  bound  for  a  for  fixed  h  =  1/N,  a  useful 
estimation  for  a  can  be  derived  as 

a  <  O.OITV  (B.IO) 

In  what  follows,  condition  (B.IO)  is  tested  by  solving  the  following  one 


197 


dimensional  convection  equation  with  and  without  added  numerical  dissipa¬ 
tion. 


ft  + fx  —  —oh* 

f{x,0)  =  (B.ll) 


The  computational  domain  is  0  <  x  <  1  with  periodic  boundary  conditions 
at  both  ends.  Although  there  are  many  choices  for  time  marching  scheme, 
we  choose  an  implicit  dual  time  stepping  scheme  with  sub-iterations.  While 
certainly  not  the  simplest,  this  scheme  is  very  representative  to  the  time 
marching  scheme  often  found  in  the  unsteady  wall  bounded  flow  simulations. 
The  derivation  of  this  scheme  is  outlined  as  follows. 


First,  the  semi-discrete  form  of  (B.ll)  with  implicit  second  order  Euler 
scheme  is 


3/n+l  _  4fn  ^  yn-l 

2At 


+ /r '  = 


(B.12) 


After  adding  a  pseudo  time  derivative  fr  to  the  left  hand  side,  (B.12)  becomes 


fr  + 


3/n+l  _  4yn  ^  yn-1 

2At 


(B.13) 


Now  we  approximate  by  first  order  Euler  scheme  fr  =  where  5/*  = 
yfe+i  _  replace  the  index  n  -I- 1  to  A:  -I- 1  in  (B.13),  the  A  form  of  the 

sub-iteration  scheme  becomes 


(—  ^ 
^At  2 At 


dx 


dx*‘ 


2At 


dx 


dx^ 


When  the  sub-iteration  converges,  i.e.  Sf*‘  =  0,  /*+^  is  taken  as  /'*+^  ,  and 
the  right  hand  side  of  (B.14)  recovers  the  (B.12).  The  pseudo  time  step  At 
can  be  properly  chosen  to  accelerate  the  sub-iteration  convergence.  In  the 
present  computation,  the  first  and  fourth  derivatives  at  the  right  hand  side  of 
(B.14)  are  approximated  using  (B.2)  and  (B.4).  The  fourth  derivative  at  the 
left  hand  side  of  (B.14)  is  approximated  by  the  following  five  point  second 
order  accurate  scheme 


d*f  1 

~  ^  ~  +  /t+i)  +  6/i].  (B.15) 

which  preserves  the  penta-diagonal  structure  of  coefficient  matrix  at  the  left 
hand  side  of  (B.14).  The  total  number  of  grid  points  is  taken  JV  =  50,  time 
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step  At  =  0.001,  and  the  pseudo  time  step  At  =  4At.  Five  sub-iterations  are 
used  for  each  time  step.  The  comparison  of  the  numerical  solutions  with  and 
without  added  dissipation  term  are  shown  in  figure  B.2  through  B.5.  Notice 
the  small  waves  at  the  trailing  edge  of  the  numerical  solution  are  the  result 
of  dispersion  error  in  the  center  differencing  scheme  —  the  traveling  speeds 
of  these  waves  are  different  from  their  theoretical  values.  The  decrease  of  the 
peak  is  also  due  to  the  dispersion  of  the  constituent  waves  in  an  Gaussian 
function.  Although  we  have  used  the  maximum  value  of  a  =  0.5  and  a  rather 
modest  resolution  N  =  50,  the  comparison  between  the  numerical  solutions 
with  and  without  the  dissipation  at  four  and  eight  flow  through  times  is 
satisfactory.  Only  slight  damping  is  observed  on  the  side  waves;  the  large 
structure  of  numerical  solution  is  essentially  unaffected  by  the  addition  of 
the  dissipation  term  at  <t  =  0.5. 

For  different  grid  resolution.  Table  1  summarizes  the  numerical  error  (in 
£-2  norm)  with  a  =  0  and  a  =  0.5  after  four  flow  through  times.  As  can  be 
seen,  at  N  <  50,  the  error  is  smaller  for  a  =  0.5  compared  to  the  no  dissi¬ 
pation  case  because  the  unsolved  scales  are  more  damped  by  the  dissipation 
term.  Beyond  N  =  50,  the  dissipation  term  shows  essentially  no  effect  on 
the  numerical  accuracy.  The  corresponding  solutions  for  N  =  100  plotted  in 
figure  B.6  shows  no  discernible  differences  between  a  —  0  and  a  =  0.5.  In 
addition,  due  to  the  increased  resolution,  the  dispersion  error  is  significantly 
reduced  so  that  the  numerical  solutions  are  in  good  agreement  with  the  exact 
solution. 

In  practical  simulations,  a  is  typically  set,  often  by  trial  and  error,  at 
the  minimum  value  that  can  stabilize  the  computation.  This  criterion,  while 
necessary,  is  not  sufficient  to  maintain  the  formal  accuracy  of  the  original 
equation.  To  do  so,  a  should  be  further  limited  below  its  upper  bound.  The 
estimation  of  the  upper  bound  of  a  also  offers  a  better  assessment  on  the  grid 
resolution:  if  the  maximum  a  is  still  insufficient  to  stabilize  the  computation, 
it  suggests  that  the  resolution  is  critically  low  and  increasing  the  number  of 
grid  points  is  essential.  Increased  resolution,  on  the  other  hand,  also  increases 
the  upper  bound  of  o  according  to  (B.IO),  so  on  a  refined  mesh,  o  can  take 
a  larger  value  and  the  numerical  method  can  achieve,  without  sacrificing  the 
accuracy,  an  improved  robustness. 
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Number  of  grid  points  N 

No  Dissipation  a  =  0 

Dissipation  at  a  =  0.5 

25 

0.1362 

0.1009 

50 

0.0214 

0.0209 

100 

0.0025 

0.0026 

200 

0.0015 

0.0015 

Table  B.l:  L2  norm  of  the  numerical  error  after  four  flow  through  times. 


Figure  B.2:  Comparison  between  numerical  and  exact  solutions  of  an  initial 
Gaussian  function,  after  4  flow  through  time  without  numerical  dissipation. 
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20] 


Figure  B.4:  Comparison  between  numerical  and  exact  solutions  of  an  initial 
Gaussian  function,  after  8  flow  through  time  without  numerical  dissipation. 
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-0.2^ 

0 


0.2  0.4 


Figure  B.5:  Same  as  above,  but  v 


2i 


Figure  B.6:  Comparison  between  numerical  and  exact  solutions 
Gaussian  function,  after  4  flow  through  time  without  numerical 
N=100. 


of  an  initial 
dissipation, 
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Appendix  C 

Compressible  Boundary  Layer 
Over  a  Leading  Edge 

C.l  Introduction 

The  self-similar  boundary  layer  momentum  and  energy  equations  in  the  re¬ 
gion  of  the  stagnation  line  of  a  swept  cylinder  in  high  speed  compressible 
flows  were  first  obtained  by  Roshotko  and  Beckwith  ?  using  the  Stewart- 
son’s  transformation  ?.  The  solutions  to  these  equations  play  the  same  role 
in  a  two  dimensional  stagnation  point  flow  as  the  Blasius  solution  in  a  flat 
plate  boundary  layer.  As  exact  solutions  to  the  Navier-Stokes  equations,  they 
have  been  used  in  the  present  study  to  compare  with  numerical  solutions  to 
determine  the  accuracy  of  the  numerical  algorithm.  The  derivation  of  self¬ 
similar  boundary  layer  equations,  somehow  unavailable  in  the  literature,  is 
provided  in  this  note  for  completeness. 

C.2  Boundary  Layer  Approximation 

Consider  compressible,  viscous,  and  heat  conducting  gas  moving  at  an  uni¬ 
form  supersonic  speed  M^o  >  1  towards  an  infinitely  long  cylinder  plaeed  at 
a  swept  angle  0,  as  sketched  in  figure  C.l.  We  define  the  coordinates  system 
such  that  the  ar-coordinate  is  the  distance  along  the  cylinder  surface  mea¬ 
sured  in  chordwise  direction  from  the  leading  edge  stagnation  line,  y  is  the 
coordinate  normal  to  the  cylinder  surface,  and  z  is  the  spanwise  direction. 
The  subscripts  oo,  1,  and  e  denote  the  flow  quantities  at  upstream  infinity. 
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TfTT  boundary  layer  in  the  region  of  the  stagnation  line 

of  a  swept  cylinder  placed  in  a  supersonic  flow. 

immediately  downstream  the  bow  shock,  and  the  edge  of  the  boundary  layer 
The  boundary  layer  thickness  is  assumed  to  be  small  in  the  stagnation  re¬ 
fill  of  curvature  of  the  surface.  Note  that 

tpr  normal  direction  is  reduced  to  subsonic  af¬ 

ter  the  bow  shock,  the  velocity  in  spanwise  direction  w,  not  affected  by  the 
shock,  may  still  be  supersonic,  i.e. 

With  the  usual  boundary  layer  approximation,  and  noting  that  all  the 
spanwise  derivatives  are  identically  zero,  i.e. 

d  _ 

Tz  =  “  (c.i) 

the  dimensional  continuity,  momentum,  energy  and  state  equations  reduce 


Continuity; 


Momentum: 


^(/^)  ,  d{fyv) 
dx  dy 

du  du  dv  8  8n 

dx  ^  dy  dx^  dy^^dy' 


(C.2) 


(C.3) 
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(C.4) 


dw  dw  d  ,  dw. 


Energy: 


dH  dH  d  ,  n  dH ,  d  ,  ,1  d  y  +  w"^ 

^  dy  ~  dy^Pr  dy^  dy^^^Pr  dy^  2 


(C.5) 


)]  (C.6) 


State: 

p  =  pBT  (C.7) 

where  p  is  density,  u,  v,  and  w  are  the  velocity  components  in  x,  y,  and  z 
directions,  p  is  pressure;  p  is  the  dynamic  viscosity,  and  Pr  is  the  Prandtl 
number 


Fr  =  i^ 

K 

The  H  in  (C.6)  is  the  total  and  stagnation  enthalpy  defined  as 
H  =  +  +  w^)  +  E  +  ^  =  +  CpT 

It  pi 


(C.8) 


(C.9) 


where  T  is  the  temperature,  E  is  the  internal  energy,  k  is  the  heat  conduction 
coefficient  and  Cp  =  is  the  specific  heat  at  constant  pressure. 

The  continuity  and  momentum  equations  follow  directly  from  the  bound¬ 
ary  layer  approximation  and  (C.l).  The  energy  equation  (C.6)  is  derived  as 
follows.  Following  ?  (p.l58),  the  general  energy  equation  is 


DH  d  ,  1 .  -  a  ,  ar, 

where  the  strain  Cjj  and  the  dilatation  e  are 


_  \  .dui  duj. 
2^dXj  ^  dxi 


A  —  Cii 


(C.IO) 


(C.ll) 


By  the  Einstein  summation  convention  and  (C.l),  (C.IO)  becomes 


d  r  ,du  1  ..  V ,dv  du.  wdw,^ 

-{2M«(^-3A)  +  5(^  +  ^)  +  -^1} 
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+ 


+ 


-3^)+ 2 

a  ,  ar,  a ,  ar, 

d  d  2  v^  ty2  du  2  ,  aT 

a^W^(“  +T  +  T*  +  '’^“3“^l  +  '‘to) 
d .  d  ,u^  2  au  2  ,  ar 

+  "  +  y)  +  "&  -  s"^!  +  ''■^> 


Apply  the  boundary  layer  approximation,  i.e.  neglecting  terms  of  the  form 

‘  ‘  ’])  ■  ■  ■]}>  high  order  terms  of  v,  we  obtain 

the  energy  equation  as  in  (C.6) 


dH  dH 


a  e  y  w\  ar. 

dy^^dy^2  2^‘^''ai,^ 


Ar 

dy  ^dy 


)+«|.(/,_!f!±J£!), 

Cpdy^  2 


—  \JL^IL^ _ ^r  rA  + 

dy^Prdy^  dy^^^Pr 


The  velocity  and  thermal  boundary  condition  are 
At  y  =  0 


At  y  -y  00 


u  =  V  =  w  =  0  and  H  =  H„ 


or 


dl^ 

dy 


=  0 


U  =  Ue,  W  =  We,  and  H  =  He 

The  viscosity  //  is  assumed  to  be  a  linear  function  of  temperature 


(C.12) 

(C.13) 


(C.14) 


where  subscript  w  denotes  the  quantities  at  the  surface.  is  a  known 
function  of  x,  and  ^  may  be  taken  as  any  desired  function  of  T,„,  such 
as  the  Sutherland  law.  Outside  the  boundary  layer,  the  chord  wise  velocity 
satisfies  the  steady  Euler  equation 


Pe«e 


dUe 

dx 


dpe 

dx 


dp 

dx 


(C.15) 


since  p  is  a  constant  in  the  direction  normal  to  the  surface  resulting  from 
p-momentum  equation. 


208 


C.3  Stewartson  Transformation 

Define  stream  function  ip  such  that 

=  (C.16) 

Po  po 

The  Stewartson  transformation  introduces  the  transformed  coordinates  X 
and  Y  by 


^  Jo  T„,  Po  Oo  Po 

=  ^  r  -T-iy 

(^0  Jo  Po 


(C.17) 


(C.18) 


here  the  subscript  0  denotes  the  total  or  stagnation  value  of  the  corresponding 
quantity,  and  a  is  the  speed  of  sound.  From  (C.17)  and  (C.18),  we  have 


^  ^dY 

9  _  vA 

dy  ”dY 

^  _  y2A-  +  y 

dy^  ^  ay  2  ^ 


.  A 

yygy 


Furthermore,  define  the  transformed  velocities  U  and  V  as 

U  =  ipY,  V  =  -V’x 

Then,  in  terms  of  U  and  V,  we  have 

u  =  — U 

ao 


u  =  -^{UY,-VX,) 

P 


(C.19) 

(C.20) 

(C.21) 


(C.22) 


(C.23) 

(C.24) 


therefore 


pup  +  ,«,^  =  P^U^(^V)-po(VY,-VX,)^(^U) 
^  dx  dy  do  dx  ao  ^o 

=  p(^fX^[UxU  +  UYV  +  —U^] 

do 


(C.25) 
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dy 


O't  P  d 

aopQdY^^aJ  poW 


®e  \3  ^^^Pe 


^  T  ^ 

Oo  Jw^ 


Uyy 


(C.26) 


where  in  the  last  expression,  (C.14)  and  (C.4)  have  been  used.  To  derive  the 
pressure  gradient,  we  first  note  that  in  external  flow.  If  is  constant;  thus 


Now  introducing  A 


7-1 

2 


(C.27) 


,Tno  l  +  ^Mlcos'^$ 

as  the  ratio  between  the  total  temperature  in  the  streamwise  direction  and 
the  total  temperature  upstream,  we  have 


Oq  “ 


7-1 


’WZ 


alX  =a^  + 


7-1 


‘UZ 


Differentiating  it  with  respect  to  x  yields, 


1  <lQe  _  7  —  1  due 

flc  dx  2al  ^^~dx 

Since  u*  =  and  14  is  only  a  function  of  X,  we  have 


(C.29) 


(C.30) 


u. 


dUe 

dx 


flo  flo  dx 


+  ^U,xX,] 
do 


=  (— )rrf  ^7  I  *^er. 


that  is 

dUe  (f)" 

dx  “  1  +  (C.31) 

Combining  (C.28).  (C.29)  and  (C.32),  the  pressure  gradient  becomes 


dp 

dx 


A  Oq 


(C.32) 
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Substituting  (C.25),  (C.26)  and  (C.32)  into  {C.3)  yields  the  x-momentum 
equation 

UUx  +  VUy  =  [  UeUex  -  —U^  ]  +  t'o  UyY  (C.33) 

p  A  Gq 


where  uo  =  /io/Po-  However, 


1  dcie  1  ^2^2 


1  (/Ug  1  /^0\2  2 

-r“  —  (— j  ^ 


2g2  ^  dx  ^ae 


.2^d(^yu,U^X.]±(^)V 
2al  A  Go  A*  G« 


2A  Gq 

and  also  noticing  that  Pe%  =  pT^  so  {C.33)  becomes 


(C.34) 


UUx  +  VUy  =[^  +  ]  T 

Ig  2ai  A  Co 


Similarly,  we  have 

dw  dw 

=  Au(wxX.  +  wyY.)  -  p[UY.  -  vx,)-^^ 
do  Po^o 

Wy 

=  p{~)^x  {Uwx  +  V  Wy) 

Go 

(C.36) 

d  .  dw. 

Oe  P  5  r  /flex  P  dw, 

Oo  Po  dY  ^  (iQ  pQ  dY 

(C.37) 

/Oe->2  PwPe 

“  ^  'T  n2 

tt/Po 

(C.38) 

which  leads  the  spanwise  momentum  equation  (C.5)  to  (C.5) 

U Wx  +  Wy  =  Uq  Wyy 

(C.39) 

Since  we  can  also  write 

pu^  +  po^  =  Au(HxX,  +  HyY.)- p(VY.-VX,) 
OX  oy  Co 
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5  1^ 


&y^PTdy’ 


5y[^(pr  2 


P{^)X^{UHx-\-VHy) 

do 


(C.40) 


_L^A_^f  /®esp5/r 

Praopo51^^^'ao'^po9>^^ 

P  (  ®e  \2  PwPe  tt 

PrW 


Pr-1 

Pr 


^  +  w^  ^  d 

dY 


t  ''  ~r  X  ^  t  IP  \  ,  -ip^  ^  ,U^  +  W^ 

I  av'  {  o  ^ ( 


2Pr  ^  {^^)yy  +  {w^)yy  ] 


2 

t2\ 


dY^ 


(C.41) 

2 

-)] 


2 

(C.42) 


iyip^y)  -  0  has  been  used.  Thus,  the  energy  equation  {C.6)  now 
takes  the  form 


UH,  +VHy=  I  ^V\y  +  (^)yy  1  +  ^Hyy 

Now,  defining  nondimensional  functions  g  and  0  as 


9 

9 


w 


We 

h-h^ 

He-H^ 


and  because  of 

we  have 


This  gives 


Hq  =  CpTo  =  Cj,Te  +  -{ul  +  wl)  =  He 


To{l-^)9  =  (To-T^) 


H  -  H. 


To'  '  ”  -“"He-H^ 

=  T  +  —  +  ^g'^-T 

2cj,  2Cp^  ^ 

=  T  +  ^  +  To{1-X)9^-T^ 


■'0  To  2  a2  j 


(C.43) 

(C.44) 

(C.45) 

(C.46) 

(C.47) 


(C.48) 


(C.49) 
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and  further 


^  (C-“) 

ie  ■^fle  -*e  -*0  ^0 

Substitute  (C.50)  into  (C.35),  the  x-momentum  equation  becomes 

UUx  +  VUy  =  l(i  -  ^)0  -  (i  -  1)9"  +  U,V,X  +  -o  Uyy  (C.51) 
Aim  A  I  m 

Substituting  the  definitions  of  A  into  (C.28)  and  g  and  6  into  (C.39),  and 
(C-43),  we  thus  obtain 
x-momentum: 

VVx  +  VUy  =  UJJex  ( 1  +  (;^  - 1)(1  -9")  +  (^  - 1)(;^)(1  -«)]  +  "o  t/yr 

J^NO  ■»  0  -I  NO 

(C.52) 

^-momentum: 

Ugx  +  V  gy  =  I'QgvY  (C.53) 

energy: 

UOx  +  V  9y  —  -^{^YY  -  C  2  )  ■^i — I  +  i^l)i9^)yy  ]} 

t'r  Oq  i  —  oo 

(C.54) 

C.4  Similar  Solutions 

To  obtain  similar  solutions,  an  external  flow  of  the  Falkner-Skan  type 


Ue  =  CAT”*  (C.55) 

is  assumed.  Introduce  the  usual  boundary  layer  transformation 

9  =  g{v)  (C.57) 

e  =  e{ri)  (C.58) 

U  =  (C.59) 
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the  momentum  and  energy  equations  in  terms  of  /,  g,  and  6  take  the  form 

r +fr= »[}'>  - 1  -  -  ,)(j  - _  ,)(i  _  tf),  (c,6o) 

7^0  7b  /j  V  / 


ft  m  » 

g  +fg'  =  0 

+  (1  - 


where 


m  + 1 

The  boundary  condition  for  the  system  equations  are 
at  77  =  0 

f=r=g=0=O 


(C.61) 

{C.62) 

{C.63) 

(C.64) 


at  77  —>  00 

f  =  g  =  0  =  l  (C.65) 

Note  the  right  hand  side  of  the  energy  equation  contains  ite  which  may  in 
general  depend  on  x.  For  similar  solution  to  be  permissible,  this  dependence 
must  be  identically  vanish,  i.e.  the  right  hand  can  only  be  a  function  of  77. 

For  gas  with  7^1,  this  can  be  realized  if  one  of  the  following  conditions  is 
satisfied. 


•  Pr  =  1 


•  «e  =  constant 

• 

In  such  cases,  the  self-similar  solutions  can  be  found  numerically.  For  the 
details  of  the  numerical  solutions,  see  ?. 
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Appendix  D 


Vorticity  boundary  conditions 


Here  we  derive  the  second  boundary  condition  (6.59)  for  the  vorticity  Umn- 
For  simplicity,  the  subscript  mn  will  be  dropped  from  the  disturbance  quanti¬ 
ties.  Hence,  the  linearized  governing  equations  for  u,  u  and  v  can  be  written 
in  a  general  form  as 

jCu;(^,  cj)  =  0,  Cu{<f>,v,u)  =  0,  jC,„{u,u},v)  =  0  (D.l) 

where  the  expressions  for  the  Cs  are  those  in  (6.15),  (6.16)  and  (6.57),  re¬ 
spectively.  Following  a  decomposition  for  the  mean  Hiemenz  velocity  (f> 

+  =  {r)-5d)  +  (t>^  (D.2) 

the  operators  £«  and  can  also  be  decomposed  into 

+  (D.3) 

where  the  superscript  p  denotes  the  operator  in  which  (j)  has  been  replaced 
by  its  potential  form  (fP,  and  b  denotes  the  complementary  operator  resulting 
from  this  decomposition;  the  effect  of  the  Hiemenz  boundary  layer  is  thus 
incorporated  in  the  operators.  The  disturbances  w,  u  and  v  can  also  be 
naturally  decomposed  as 

oj  =  u  =  +  u^;  v  =  if  +  v'’  (D.4) 

For  the  p  quantities,  the  governing  equations  are 

Ci{r,uf)  =  0-,  Cl{r,v^,un  =  0-,  C„{uruf,if)  =  0  (D.5) 
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The  governing  equations  for  the  b  quantities  follow  directly  from  the  decom¬ 
positions  in  (D.3)  and  (D.4).  For  velocities  u’*  and  v^,  we  enforce  the  same 
boundary  conditions  as  those  for  the  original  u  and  v,  i.e. 


vF  =  0  at  Tf  =  0  and  77  ->  oo 


t;’’  =  0  and 


dv^ 

-7—  =  0  at  n  =  0 
dr) 


For  the  tu*’,  (6.56)  leads  to  the  first  boundary  condition 


(D.6) 

(D.7) 


^  U)oo  as  T)  00 


{D.8) 


where  cJqo  is  the  initial  disturbance  vorticity  introduced  far  upstream.  The 
second  boundary  condition  for  cJ'  can  be  derived  from  the  following  fact  that 
asT)  00 

„B  ,  /  1  for  fundamental  mode 

"“"jo  else  (D-9) 

which  is  implied  by  (6.21).  To  see  this,  solving  the  last  equation  in  (D.5) 
subject  to  the  boundary  condition  (D.7),  we  obtain 


gn*o»?  p  p-nkon  rrj 

^  = - Y-  +  dr]'  +  I  {iujP  drf 

(D.IO) 

For  (D.9)  to  be  realizable,  must  remain  bounded  as  r)  00.  So  the 
coefficient  of  in  (D.IO)  must  go  to  zero  as  77  — >  00,  i.e. 


prj 


poo 

/  {iuf  +  df)’  =  0 

^0 


Moreover,  notice  that 


uP  =  0 


(D.ll) 


(D.12) 


as  a  result  of  (<^)  =  0  in  and  the  homogeneous  boundary  condition 
(D.6).  Thus  (D.ll)  reduces  to 


poo 

I  ^-nkor/  _  Q 

Jo 


(D.13) 


which  serves  as  the  second  boundary  condition  for  cjP.  Since  the  general 
expression  for  cjP  has  been  obtained  in  (6.46),  (D.8)  and  (D.13)  can  thus  be 
used  to  specify  the  two  arbitrary  constants  therein. 
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Once  ct/P  and  v’*  are  known,  the  corresponding  b  quantities  can  be  readily 
solved.  Notice  that,  by  construction,  v!’  and  satisfy  homogeneous  boundary 
conditions,  and  the  boundary  condition  for  u}^  for  7;  — >  oo  is  also  homoge¬ 
neous. 
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Appendix  E 

Vorticity  asymptotes 


The  general  asymptotic  expression  for  the  confluent  hyper-geometric  function 
with  large  real  argument  is  (Abramowitz  and  Stegun,  1970) 

^(a;  c;  -z)  =  [1  -t-  0(1^1“^)]  as  ^ oo;  ph(z)  =  0  (E.l) 

Since  the  inflow  boundary  is  assumed  to  be  far  upstream,  i.e.  Hq  1,  the 
asymptotic  expressions  for  Mi  and  M2  in  (6.60)  are  as  follows. 


Ml 


M2 


_ ^^^2^  (Ih.)^-l+iTn<7o 


(E.2a) 


(E.2b) 


When  A  — >  00,  by  the  nature  of  F  function  on  complex  plane,  |r(z)|  decrease 
rapidly  along  the  imaginary  axis.  So  (E.2)  shows  that  for  large  A,  M2  > 
Ml.  Moreover,  |/i|  and  |/2(  in  (6.62)  can  be  shown  of  the  same  order  of 
magnitude,  thus  compared  to  \I1\I\I2\M2,  Mi  can  be  neglected.  Substituting 
these  relations  into  (6.61)  yields  the  asymptotic  expression  for  the  amplitudes 
of  Cmn  and  at  large  A: 


D„ 


^/2^(5)  h  'v/2' 

1(2  •*•  «  +  ■?)  -p  3. , 


218 


In  order  to  evaluate  Wmn  at  the  wall,  the  amplitudes  of  Cmn  and  Dmn  may 
be  estimated  more  explicitly  if  the  fundamental  frequency  gq  is  low.  In  fact, 
for  the  r  function  of  a  complex  argument,  it  follows 


(E.4) 


and  as  <To  1,  it  becomes 


^,1  V?  .mtro. 

"^(2  +  21+'— ' 


r(i  +  i) 

\/l  + 


(E.5) 


where 


Substitute  these  expressions  into  (E.3a),  and  expand  the  F  function  in  terms 
of  power  series  of  ^  and  rn^ol  up  to  the  first  order,  we  have 


\Crr 


IMXi+I) 


\h\^/2T{\)yJl  +  a^  ^ 


n^,r.  Q!„m"  2i 


~  /fo[l  +  (ai-lnFo)y][l-^a2] 


(E.6) 


where  ci  =  |(ln2  +  7)  and  7  =  0.5772156--  -  is  the  Euler  constant.  The 
expression  for  \Dmn\  can  be  similarly  obtained.  Notice  that  the  value  of 
I/1I/I/2I,  as  mentioned  before,  is  slow-varying  and  of  order  one.  For  instance, 
in  the  limit  of  low  frequency  and  large  scale,  i.e.  oo  — >  0  and  A  — >  00  , 

h  ~  y/fll-erf(^));  h  ~  ^Ei(l.i)  (E.7) 

where  erf  and  Ei  are  the  error  function  and  exponential  integral,  we  have 
therefore  |/2|/|/i|  ->  0.70378178- • To  simplify  the  discussion,  the  depen¬ 
dence  of  /1//2  on  A  and  gq  will  be  neglected,  i.e.  |/i|/|72|  is  treated  as  a 
constant. 


219 


Bibliography 


R.  Abid  and  C.  G.  Speziale.  The  freestream  matching  condition  for  stagna¬ 
tion  point  turbulent  flows  -  an  alternative  formulation.  Technical  Report 
AM-94-001,  Boston  University,  1994. 

M.  Abramowitz  and  I.  A.  Stegun,  editors.  Handbook  of  Mathematical  Func¬ 
tions  with  Formulas,  Graphs,  and  Mathematical  Tables.  Dover  Publica¬ 
tions,  Inc.,  New  York,  1970. 

F.  E.  Ames.  The  influence  of  large-scale  high-intensity  turbulence  on  vane 
heat  transfer.  ASME  Journal  of  Turbomachinary,  199:23-30,  1997. 

F.  E.  Ames  and  R.  J.  Moffat.  Heat  transfer  with  high  intensity,  large  scale 
turbulence:  The  flat  plate  turbulent  boundary  layer  and  the  cylindrical 
stagnation  point.  In  Dept,  of  Mech.  Engng.  Report  No.  HMT-44,  Stanford 
University,  Stanford,  CA,  1990. 

F.  E.  Ames,  C.  Wang,  and  P.  A.  Barbot.  Measurement  and  prediction  of 
t  e  influence  of  catalytic  and  dry  low  NOx  combustor  turbulence  on  vane 
surface  heat  transfer.  In  Proceedings  of  ASME  TURBO  EXPO  2002,  GT- 
2002-30524,  Amsterdam,  The  Netherlands,  June  3-6  2002. 

B.  Andreotti,  S.  Douady,  and  Y.  Couder.  An  experiment  on  two  aspects  of 

the  interaction  between  strain  and  vorticity.  J.  Fluid  Mech.,  444:151-174, 
2001  • 

S.  Bac,  S.  K.  Lele,  and  H.  J.  Sung.  Influence  of  inflow  disturbances  on 
stagnation-region  heat  transfer.  ASME  Journal  of  Heat  Transfer,  122: 
258—265,  2000. 


220 


S.  Bae,  S.  K.  Lele,  and  H.  J.  Sung.  Direct  numerical  simulation  of  stagna¬ 
tion  region  flow  and  heat  transfer  with  free-stream  turbulence.  Physics  of 
Fluids,  15(6);1462-84,  2003. 

J.E.  Bardina,  P.  G.  Huang,  and  T.  J.  Coakley.  'Birbulence  modeling  valida¬ 
tion.  AIAA  Paper  97-2121, 1997. 

G.  K.  Batchelor.  An  introduction  to  fluid  dynamics.  Cambridge  university 
press,  1967. 

G.  K.  Batchelor  and  I.  Proudman.  The  effects  of  rapid  distortion  of  a  fluid 
in  turbulent  motion.  Q.  J.  Mech.  Appl.  Maths,  7:83-103,  1954. 

R.  M.  Beam  and  R.F.  Warming.  An  implicit  factored  scheme  for  the  com¬ 
pressible  Navier-Stokes  equations.  AIAA  Journal,  16(4):393-402,  April 
1978. 

M.  Behnia,  S.  Parneix,  Y.  Shabany,  and  P.  A.  Durbin.  Numerical  study  of 
turbulent  heat  tranfer  in  confined  and  unconfined  impinging  jets.  Int.  J. 
Heat  Fluid  Flow,  20:1-9, 1999. 

J.  Botcher  and  E.  Wedemeyer.  The  flow  downstream  of  screen  and  its  in¬ 
fluence  on  the  flow  in  the  stagnation  region  of  cylindrical  bodies.  J.  Fluid 
Mech.,  204:501-522,  1989. 

W.  R.  Briley  and  H.  McDonald.  Solution  of  the  three-dimensional  compress¬ 
ible  Navier-Stokes  equations  by  an  implicit  technique.  In  Proceedings  of  the 
Fourth  International  Conference  on  Numerical  Methods  in  Fluid  Dynam¬ 
ics,  Lecture  Notes  in  Physics,  volume  35,  pages  105-110.  Springer- Verlag, 
Berlin,  1975. 

M.  Champion  and  P.  A.  Libby.  Asymptotic  analysis  of  stagnating  turbulent 
flows.  AIAA  J.,  29{l):16-24, 1991. 

S.  Scott  Collis.  A  Computaitonal  Investigation  of  Receptivity  in  High-Speed 
Flow  Near  A  Swept  Leading-Edge.  PhD  thesis,  Stanford  University,  April 
1997. 

G.  Comte-Bellot  and  S.  Corrsin.  Simple  Eulerian  time  correlation  of  full  and 
narrow-band  velocity  signals  in  grid  generated  ‘isotropic’  turbulences.eddy 
shocklets  in  decaying  compressible  turbulence.  J.  Fluid  Mech.,  48:273-337, 
1971. 


221 


T.  J.  Craft,  L.  J.  W.  Graham,  and  B.  E.  Launder.  Impinging  jet  sutdies  for 
turbulence  model  assesment-ii.  an  examination  of  the  performance  of  four 
turbulence  models.  Int.  J.  Heat  Mass  TVansf.,  36(10):2685-2697, 1993. 

M.  R.  Dhanak  and  J.  T.  Stuart.  Distortion  of  the  stagnation-point  flow  due 
to  cross-stream  vorticity  in  the  external  flow.  Phil.  Trans.  R  Soc  Land 
A,  352:443-452,  1995. 

K.  Dullenkopf  and  R.  E.  Mayle.  An  account  of  free-stream-turbulence  length 
scale  on  laminar  heat  transfer.  ASMS  Journal  of  Turbomachinaru  117- 
401-406,  1995. 

P.  Durbin.  On  the  A:  -  c  stagnation  point  anomaly.  Int.  J.  Heat  and  Fluid 
Flow,  17:89-90, 1996. 

C.  A.  J.  Fletcher.  Computational  techniques  for  fluid  dynamics.  Berlin  ;  New 
York  :  Springer- Verlag,  2nd  edition,  1991. 

N.  Frossling.  Evaporation  heat  transfer  and  velocity  distribution  in  two- 
dimensional  and  rotationally  symmetrical  laminar  boundary-layer  flow. 
Technical  Report  TM-1432,  NACA,  1940. 

M.  Germano.  Averaging  invariance  of  the  turbulent  equations  and  similar 
subgrid-scale  modeling.  Manuscript,  116,  Center  for  Turbulence  Research, 


M.  M.  Gibson  and  B.  E.  Launder.  Ground  effects  on  pressure  fluctuation  in 
the  atmospheric  boundary  layer.  J.  Fluid  Mech.,  86:491-591,  1978. 

W.  H.  Giedt.  Investigation  of  variation  of  point  unit  heat  transfer  coefficient 
around  a  cylinder  normal  to  an  airstream.  Trans.  ASME  71  pages  375- 
381,  1949. 

M.  E.  Goldstein.  Unsteady  vortical  and  entropic  distortions  of  potential  flows 
round  arbitrary  obstacles.  J.  Fluid  Mech.,  89:433-468,  1978. 

R.  J.  Goldstein,  editor.  Heat  transfer  in  gas  turbine  systems,  volume  934. 
Annals  of  the  New  York  Academy  of  Scienes,  2001. 

H.  Gortler.  Three  dimensional  instability  of  the  plane  stagnation  flow  with 
respect  to  vortical  disturbance.  In  H.  Gortler  and  W.  Tollmien,  editors. 


222 


Fifty  years  of  boundary  layer  research,  pages  304-314.  Vieweg  and  Sohn, 
1955.  Braunschweig. 

C.  E.  Grosch  and  H.  Salwen.  Oscillating  stagnation  point  flow.  Proc.  R.  Soc. 
Land.  A,  384:175-190,  1982. 

G.  Hammerlin.  On  instability  theory  of  plane  stagnation  point  flow.  In 
H.  Gortler  and  W.  Tollmien,  editors,  Fifty  years  of  boundary  layer  research, 
pages  315-327.  Vieweg  and  Sohn,  1955.  Braunschweig. 

P.  E.  Hancock  and  P.  Bradshaw.  The  effect  of  free-stream  turbulence  on 
turbulent  boundary  layers.  ASME  Journal  of  Fluids  Engineering,  105: 
284-289,  1983. 

B.  G.  Van  Der  Hegge-Zijnen.  Heat  transfer  from  horizontal  cylinders  to  a 
turbulent  air  flow.  Appl.  Sci.  Res.  A,  7:205-223,  1957. 

J.  C.  R.  Hunt.  A  theory  of  turbulent  flow  round  two-dimensional  bluff  bodies. 
J.  Fluid  Mech.,  61:625-706,  1973. 

J.  C.  R.  Hunt  and  J.  M.  R.  Graham.  Free-stream  turbulence  near  plane 
boundaries.  J.  Fluid  Mech.,  84:205-235,  1978. 

Y.  H.  Im,  K.  Y.  Huh,  and  K-Y.  Kim.  Analysis  of  impinging  and  coun¬ 
tercurrent  stagnating  flows  by  Reynolds  stress  model.  ASME  J.  Fluid 
Engineering,  124:706-718,  2003. 

E.  Isaacson  and  H.  B.  Keller.  Analysis  of  numerical  methods.  Dover  publi¬ 
cation,  Inc.,  1993. 

M.  Ishigaki.  Periodic  boundary  layer  near  a  two-dimensional  stagnation 
point.  J.  Fluid  Mech.,  43:477-486, 1970. 

M.  Kato  and  B.  E.  Launder.  Modelling  flow-induced  oscillation  in  turbulent 
flow  around  a  square  cylinder.  ASME  FED,  157:189-199, 1993. 

O.  S.  Kerr  and  J.  W.  Dold.  Periodic  steady  vortices  in  a  stagnation-point 
flow.  J.  Fluid  Mech.,  276:307-325, 1994. 

J.  Kestin.  The  effect  of  free-stream  turbulence  on  heat  transfer  rates.  In  T.  F. 
Irvine  and  J.  P.  Hartnett,  editors.  Advances  in  heat  transfer,  volume  3. 
Academic  Press,  1966. 


223 


J.  Kestin,  P.  F.  Madder,  and  H.  H.  Sogin.  The  influence  of  turbulence  on  the 
trnafer  of  heat  to  cylinders  near  the  stagnation  point.  ZAMP,  12:115-132 
1961. 

J.  Kestin  and  R.  T.  Wood.  On  the  stability  of  two-dimensional  stagnation 
flow.  J.  Fluid  Mech.,  44:461-479,  1970. 

J.  Kestin  and  R.  T.  Wood.  The  influence  of  turbulence  on  mass  transfer  from 
cylinders.  ASMS  J.  Heat  Transfer,  930:321-327, 1971. 

S.  Lee,  S.  K.  Lele,  and  P.  Moin.  Eddy  shocklets  in  decaying  compressible 
turbulence.  Physics  of  Fluids,  A  3(4):657-664,  1991. 

S.  Lee,  S.  K.  Lele,  and  P.  Moin.  Direct  numerical  simulation  of  isotropic 
turbulence  interacting  with  a  weak  shock  wave.  J.  Fluid  Mech  251-533- 
562,  1993. 

S.  K.  Lele.  Compact  finite  difference  schemes  with  spectral-like  resolution. 
Journal  of  Computational  Physics,  103(l):16-42,  1992. 

M.  J.  Lighthill.  The  response  of  laminar  skin  friction  and  heat  transfer  to 
fluctuations  in  the  stream  velocity.  Proc.  Roy.  Soc.  A,  224:1-23,  1954. 

N.  Lin.  Receptivity  of  the  boundary  layer  over  a  fiat  plate  xvith  different 
leading-edge  geometries:  Numerical  simulations.  PhD  thesis,  Arizona  State 
University,  1992. 

R.  S.  Lin  and  M.  R.  Malik.  On  the  stability  of  attachment-line  boundary 
layers.  Part  1.  the  incompressible  swept  Himenz  flow.  J.  Fluid  Mech  31T 
23^255,  1996. 

G.  W.  Lowery  and  R.  I.  Vachon.  Effect  of  turbulence  on  heat  transfer  from 
heated  cylinders.  Int.  J.  Heat  Mass  Transfer,  18(No.  11):1229-1242, 1975. 

C.  Lui.  A  numerical  investigation  of  shock  associated  noise.  PhD  thesis, 
Stanford  University,  2003. 

C.  Lui  and  S.  K.  Lele.  Direct  numerical  simulation  of  spatially  developing 
compressible  turbulence  mixing  layers.  AIAA  Paper  2001-0291,  2001. 

M.  J.  Lyell  and  P.  Huerre.  Linear  and  nonlinear  stability  of  plane  stagnation 
flow.  J.  Fluid  Mech.,  161:295-312,  1985. 


224 


G.  Medic  and  P.  A.  Durbin.  Toward  improved  prediction  of  heat  transfer  on 
turbine  blades.  ASMS  J.  of  Turbomachinery,  124:187-192,  2002. 

A.  B.  Mehendale,  J.  C.  Han,  and  S.  Ou.  Influence  of  high  mainstream  tur¬ 
bulence  on  leading  edge  heat  transfer.  ASME  Journal  of  Heat  Transfer, 
113:843-850, 1991. 

G.  J.  Merchant  and  S.  H.  Davis.  Modulated  stagnation-point  flow  and  steady 
streaming.  J.  Fluid  Mech.,  198:543—555, 1989. 

L.  M.  Milne-Thomson.  Theoretical  hydrodynamics.  The  Macmilian  Company, 
New  York,  4th  edition,  1960. 

P.  Moin,  K.  Squires,  W.  Cabot,  and  S.Lee.  A  dynamic  subgrid-scal  model 
for  compressible  turbulence  and  scalar  transport.  Physics  of  Fluids,  A3 
(11),  November  1991. 

M.  V.  Morkovin.  Bypass-transition  research:  issues  and  philosophy.  Tech¬ 
nical  Report  Tech.  Rep.  AFFDL-TR-68-149,  Air  Force  Flight  Dynamics 
laboratory,  Wright-Paterson  Air  Force  Base,  1969. 

M.  V.  Morkovin.  On  the  question  of  instabilities  upstream  of  cylindrical 
bodies.  NASA  Contractor  Report,  3231,  1979. 

H.  M.  Nagib  and  P.  R.  Hodson.  Vortices  induced  in  a  stagnation  region 
by  wakes.  In  L.  S.  Fletcher,  editor.  Aerodynamic  heating  and  thermal 
protection  systems,  volume  59  of  Progress  in  astronautics  and  aeronautics, 
1978. 

W.  Nakayama.  Heat-transfer  engineering  in  systems  integration:  outlook 
for  closer  coupling  of  thermal  and  electrical  designs  of  computers.  IEEE 
Transactions  on  Components,  Packaging,  and  Manufacturing  Technology, 
PaH  A,  18(4):818-26,  1995. 

A.  N.  Oo  and  C.  Y.  Ching.  Effect  of  turbulence  with  different  vortical  struc¬ 
tures  on  stagnation  region  heat  transfer.  ASME  Journal  of  Heat  Transfer, 
123:665-674,  2001. 

S.  Pameix,  M.  Behnia,  and  P.  A.  Durbin.  Predictions  of  turbulent  heat 
transfer  in  an  axisymmetric  jet  impinging  on  a  heated  pedestal.  ASME  J. 
Heat  Tranfer,  121:43-49, 1999. 


225 


T.  J.  Pedley.  Two-dimensional  boundary  layers  in  a  free  stream  which  oscil¬ 
lates  without  reversing.  J.  Fluid  Mech.,  55:359-383,  1972. 

B.  Perot  and  P.  Moin.  Shear-free  turbulent  boundary  layers,  part  1.  physical 
insights.  J.  Fluid  Mech.,  295:199-227,  1995. 

N.  A  V.  Piercy  and  E.  G.  Richardson.  The  variation  of  velocity  amplitude 
colse  to  the  surface  of  a  cylinder  moving  through  a  viscous  fluid.  Phil 
Mag.  Series  7,  6:970-977,  1928. 

N.  A.  V.  Piercy  and  E.  G.  Richardson.  The  turbulence  in  front  of  a  body 
moving  through  a  viscous  fluid.  Phil.  Mag.  Series  7,  9:1038-1041,  1930. 

accuracy  and  the  use  of  implicit  methods.  AIAA  paper 


M.  M.  Rai.  Navier-Stokes  simulations  of  blade-vortex  interaction  using  high- 
order  accurate  upwind  schemes.  AIAA  paper  87-0543,  1987. 

E.  Reshotko  and  I.  E.  Beckwith.  Compressible  laminar  boundary  layer  over 
a  yawed  infinite  cylinder  with  heat  transfer  and  arbitrary  prandtl  number 
Report  1379,  NACA,  1957. 


^  i  J-  VanFossen.  Increased  heat  transfer  to  elliptical  leading 

edges  due  to  spanwise  variations  in  the  free-stream  momentum:  Numerical 
and  experimental  results.  AIAA  paper  92-3070,  1992. 


W.  J.  Sadeh  and  H.  J.  Brauer.  A  visual  investigation  of  turbulence  in  sta- 
gantion  flow  about  a  circular  cylinder.  J.  Fluid  Mech.,  99:53-64,  1980. 

M.  Smith  and  A.  M.  Kuethe.  Effects  of  turbulence  on  laminar  skin  friction 
and  heat  transfer.  Phys.  Fluids,  9(No.  12):2337-2344, 1966. 


^  numerical  study  of  leading-edge  contamination. 

AGARDC,  438(5.1-5.13):457-468,  1989. 


P.  R.  Spalart  and  S.  R.  Allmaras.  A  one-equation  turbulence  model  for 
aerodynamic  flows.  AIAA  paper  92-0439,  1992. 

C.  G.  Speziale,  S.  Sarkar,  and  T.  B.  Gatski.  Modelling  the  pressure-strain 
correlation  of  turbulence:  an  invariant  dyanmic  systems  approach.  J.  Fluid 
Mech.,  227:245-272,  1991. 


226 


S.  P.  Sutera.  Vorticity  amplification  in  staganation-point  flow  and  its  effect 
on  heat  transfer.  J.  Fluid  Mech.,  21:513-534,  1965. 

D.  B.  Taulbee  and  L.  Tran.  Stagnation  streamline  turbulence.  AIAA  J.,  26 
(8):1011-1013, 1988. 

V.  Theofilis,  A.  Fedorov,  D.  Obrist,  and  U.  Ch.  Dallmann.  The  extended 
Gortler-Hammerlin  model  for  linear  instability  of  three-dimensional  incom¬ 
pressible  swept  attachment-line  boundary  layer  flow.  J.  Fluid  Mech..,  487. 
271-313,  2003. 

K.  A.  Thole,  R.  W.  Radomsky,  M.  B.  Kang,  and  A  Kohli.  Elevated  freestream 
turbulence  effects  on  heat  transfer  for  a  gas  turbine  vane.  Int.  J.  Heat  and 
Fluid  Flow,  23:137-147,  2002. 

R.  M.  TVaci  and  D.  C.  Wilcox.  Freestream  turbulence  effects  on  stagnation 
point  heat  transfer.  AIAA  J.,  13(7):890-896, 1975. 

G.  I.  Tsourakis,  D.  G.  Koubogiannis,  and  K.  C.  Giannakoglou.  Transition 
and  heat  transfer  predictions  in  a  turbine  cascade  at  various  free-stream 
turbulence  intensities  through  a  one-equation  turbulence  model.  Int.  J. 
Numer.  Meth.  Fluids,  38:1091-1110,  2002. 

M.  Van  Dyke.  Perturbation  methods  in  fluid  mechanics.  The  parabolic  press, 
Stanford,  California,  1975. 

G.  J.  Van  Fossen,  R.  J.  Simoneau,  and  C.  Y.  Ching.  Influence  of  turbulence 
parameters,  Reynolds  number  and  body  shape  on  stagnation  region  heat 
transfer.  ASME  Journal  of  Heat  Transfer,  117:593—603, 1995. 

G.  J.  VanFossen  and  R.  S.  Bunker.  Augmentation  of  stagnation  region  heat 
transfer  due  to  turbulence  from  a  DLN  can  combuster.  NASA/TM-200- 

210241,  2000. 

S.  Venkateswaran,  P.  E.  O.  Buelow,  and  C.  L.  Merkle.  Development  of 
linearized  preconditioning  methods  for  enhancing  robustness  and  efficiency 
of  euler  and  navier-stokes  computations.  AIAA  Paper  97-2030, 1997. 

C.  R.  Wang  and  F.  C.  Yeh.  Application  of  turbulence  modeling  to  predict 
surface  heat  transfer  in  stagnation  flow  region  of  circular  cylinder.  NASA, 
Technical  Paper,  2758,  1987. 


227 


D.  C.  Wilcox.  IVirbulence  modeling:  An  overview.  AIAA  Paper  2001-0721 
2001. 

S.  D.  R.  Wilson  and  I.  Gladwell.  The  stability  of  a  two-dimensional  stagnation 
flow  to  three-dimensional  disturbances.  J.  Fluid  Mech.,  84:517-527, 1978. 

R.  W.  Wlezien.  Measurement  of  acoustic  receptivity.  25th  AIAA  Fluid 
Dynamics  Conference,  AIAA  94-2221, 1994. 

Z.  Xiong.  Stagnation  point  flow  and  heat  transfer  under  free-stream  turbu¬ 
lence.  PhD  thesis,  Stanford  University,  2004. 

Z.  Xiong  and  S.  K.  Lele.  Numerical  study  of  leading-edge  heat  transfer  under 
free-stream  turbulence.  AIAA  Paper  2001-1016,  January  2001. 

N.  R.  Yardi  and  S.  P.  Sukhatme.  Effect  of  turbulence  intensity  and  integral 
length  scale  of  a  turbulent  free  stream  on  forced  convection  heat  transfer 
from  a  circular  cylinder  in  cross-flow.  In  Proceedings  6th  international 
heat  transfer  conference,  number  5  in  FC(b)-29,  pages  347-352,  Toronto, 
Canada,  1978. 

F.  C.  Yeh,  S.  A.  Hippensteele,  G.  J.  Van  Fossen,  P.  E.  Poinsatte,  and 
A.  Ameri.  High  reynolds  number  and  turbulence  effect  on  aerodynam¬ 
ics  and  heat  transfer  in  a  turbine  cascade.  AIAA-93-2252,  1993. 


228 


